Clear the environment

rm(list = ls())
knitr::opts_chunk$set(warnings=FALSE) 

Load packages

Before loading the packages install them using install()

library(tidyverse)
library(summarytools)
library(dplyr)
library(readr)
library(vegan)
library(patchwork)
library(pander)
library(ggeffects)
library(performance)
library(knitr)
library(car)
library(corrplot)
library(MASS)
library(MuMIn)
library(geosphere)
library(mgcv)
library(glmmTMB)
library(lme4)
library(factoextra)
library(fitdistrplus)
#If Issues installing packages
options("install.lock"=FALSE)

Part 1: Importing and examining data

In this part the data frames were inspected in order to check the variable names and types. In some situations the data set was slightly modified for further manipulation of them.

1a) Import data

Environmental1 <- read.csv("Data/Environmental_variables.csv")
Metadata <- read.csv("Data/Metadata_lat2018-2019.csv")
Herbivory <- read.csv("Data/Herbivory_latitudinal2018-2019.csv")
Benthic <- read.csv("Data/Percent_cover.csv")
Phlorotannins <- read.csv("Data/Phlorotannins.csv")
Predation <- read.csv("Data/Predation_latitudinal2018-2019.csv")
Size <- read.csv("Data/Size_estimations.csv")
Transects <- read.csv("Data/Transect2018-2019.csv")
Sal_data <-read.csv("Data/extracted_sal_data.csv")
temp_data <- read.csv("Data/extracted_temp_data.csv")
o2_data <- read.csv("Data/extracted_o2_data.csv")

1b) Inspect structure of data frames

  1. Site metadata
view(dfSummary(Metadata))

#Rename temperature variables to avoid confusion with replicate measures in Environmental data
Metadata <- rename(Metadata, Temperature_site = Temperature_C)
  1. Environmental data
view(dfSummary(Environmental1))

Environmental <- Environmental1 %>%
  rename(D_oxygen_percent = X., D_oxygen_mgL = D_oxygen) %>% #Rename oxygen percent variable
  dplyr::select(-Date,-Campaign,-Latitude) #Remove variables that are in Metadata
  1. Benthic cover
view(dfSummary(Benthic))
Benthic <- Benthic %>%
  dplyr::select(-Date) #Remove variables that are in Metadata
  1. Herbivory
view(dfSummary(Herbivory))
Herbivory <- Herbivory %>%
  dplyr::select(-Date,-Campaign,-Season) %>%# Remove variables that are in Metadata
  mutate(Consumed_proportion = (Initial_length_mm-Final_Length_mm)/Initial_length_mm)

It was necessary to create major categories for Bait

Herbivory <- mutate(Herbivory, Bait2 = case_when(Bait == "Lessonia_spicata" ~ "Lessonia sp.",
                                                 Bait == "Lessonia_berteroana" ~ "Lessonia sp."))
  1. Predation
view(dfSummary(Predation))
Predation <- Predation %>%
  dplyr::select(-Date,-Campaign,-Season)  %>%# Remove variables that are in Metadata
mutate(Consumed_proportion = (Initial_Individuals-Remaining_individuals)/Initial_Individuals, Consumed_individuals = Initial_Individuals-Remaining_individuals) #Add Consumed_proportion column

Create major categories for Bait (Squidpops vs Crabs) and Assay time (2 hours and 24 hours)

Predation <- mutate(Predation, Bait2 = case_when(Bait == "Squidpops" ~ "Squidpops",
                                                 Bait != "Squidpops" ~ "Crabs"), Assay_time2 = case_when(Assay_time == "24" ~ "24",
                                                 Bait != "24" ~ "2"))
  1. Phlorotannins
view(dfSummary(Phlorotannins))
Phlorotannins <- Phlorotannins %>%
  dplyr::select(-Date,-Latitude) # Remove variables that are in metadata 
  1. Size
view(dfSummary(Size))
Size <- Size %>%
  dplyr::select(-Date,-Campaign) # Remove variables that are in Metadata
  1. Transects
view(dfSummary(Transects))
Transects <- Transects %>%
  dplyr::select(-Date,-Campaign) # Remove variables that are in Metadata

Part 2: Manipulating datasets

In this section data sets are manipulated for analyses. New data sets with means and combined with metadata and consumption are created.

2a) Environmental data

  1. Join Environmental data with metadata
Environmental_metadata <- Environmental %>%
    left_join(Metadata, by = "Site")
  1. Join Environmental data with consumption
  • First summarise environmental data per site
Environmental_2 <- Environmental_metadata %>%
  group_by(Site, Latitude) %>%
  dplyr::select(-D_oxygen_percent, -Conductivity_Ms, -C_f, -C_i, -H_Visibility) #Keep only variables that will be used in the analysis

Environmental_2$Day_length <- daylength(Environmental_2$Latitude, Environmental_2$Date)


Environmental_sum_long <- Environmental_2 %>%
  dplyr::select(Site, Latitude, Temperature, D_oxygen_mgL, Salinity_ppt, V_Visibility, Day_length) %>%
  gather(Env_Variable, Value, Temperature:Day_length) %>%
  group_by(Site, Latitude, Env_Variable)%>%
  summarise(Mean_value = mean(Value, na.rm=T),n = n(), sd = sd(Value,na.rm=T), se = sd/sqrt(n)) 

Environmental_sum <- Environmental_sum_long %>%
  dplyr::select(Site, Env_Variable, Mean_value) %>%
pivot_wider(names_from = Env_Variable, values_from = Mean_value,
              values_fill = 0)
  • Then, add environmental data to consumption

  • Predation vs environmental data sets

Predation_env <- Predation %>%
    left_join(Environmental_sum, by = "Site")
  • Herbivory vs environmental data sets
Herbivory_env <- Herbivory %>%
  left_join(Environmental_sum, by = "Site")

2b): Benthic data

  1. Join benthic data with metadata
Benthic_metadata <- Benthic  %>%
    left_join(Metadata, by = "Site")

Benthic_metadata_wide <- Benthic_metadata %>%
  group_by(Site, Latitude, Type) %>%
  pivot_wider(names_from = Type, values_from = Coverage,
              values_fill = 0)
  1. Calculate mean, sd, and se by group
Benthic_mean <- Benthic_metadata %>%
  group_by(Site, Latitude, Type)%>%
  summarise(Mean_cover = mean(Coverage),n = n(), sd = sd(Coverage), se = sd/sqrt(n)) 
  1. Join benthic cover data with consumption
  • First, summarise benthic cover per site and change it to the wide format
Benthic_sum_wide <- Benthic %>%
  group_by(Site, Type) %>%
  summarise(cover_mean = mean(Coverage)) %>%
  pivot_wider(names_from = Type, values_from = cover_mean,
              values_fill = 0)
Benthic_sum_wide_var <- dplyr::select(Benthic_sum_wide, Site, Kelp, Bare_rock, Understorey_algae)
  • Then, add benthic data to consumption

  • Predation_env + Benthic summarised data

Predation_2 <- Predation_env %>%
    left_join(Benthic_sum_wide_var, by = "Site")
  • Herbivory_env + Benthic summarised data
Herbivory_2 <- Herbivory_env %>%
  left_join(Benthic_sum_wide_var, by = "Site")

2c): Fish length

  1. Join Size data with metadata
Size_metadata <- Size %>%
left_join(Metadata, by = c("Site"))
  1. Calculate mean, sd and se per group
Size_mean_sp<- Size_metadata %>%
  group_by(Site, Latitude, Species, Type)%>%
  summarise(mean_Size_sp = mean(Size_cm))
Size_mean <- Size_mean_sp %>%
   group_by(Site, Latitude, Type)%>%
  summarise(mean_Size = mean(mean_Size_sp),n = n(), sd = sd(mean_Size_sp), se = sd/sqrt(n)) 
  1. Calculate mean of fish size and add it to consumption
Size_mean_wide <- Size_mean_sp %>%
  filter(Type == "Fish") %>%
  group_by(Site, Type)%>%
  summarise(mean_Size = mean(mean_Size_sp)) %>%
  pivot_wider(names_from = Type, values_from = mean_Size,
              values_fill = 0) %>%
  rename(Fish_size = Fish)
  • Add fish length data to consumption

  • Predation_env + fish size data

Predation_3 <- Predation_2 %>%
    left_join(Size_mean_wide, by = "Site")

2d) Species list (abundances and richness)

  1. Add metadata to Transects data set
#Join transects with metadata
Transect_metadata <- Transects %>%
left_join(Metadata, by = c("Site"))

#Summary of species from transects
Transect_sum <- Transect_metadata %>%
  filter(Method == "Transect") %>% #Select only data taken with visual transect method
  group_by(Site, Latitude, Census, Species) %>%
  summarise(Total_abundance = sum(Individuals))
  1. Calculate abundances
  • Calculate abundances by species

First, change data frame to the wide format

Transect_wide <- Transect_sum %>%
   ungroup()  %>%
  pivot_wider(names_from = Species, values_from = Total_abundance,
              values_fill = 0)

Convert back to long format with Species and Abundance as columns. This will add rows with zeroes for all the species missing in each video

Transect_long <- Transect_wide %>%
  gather(Species, Abundance, Anisotremus_scapularis:Labrisomus_philippii)

#Obtain mean abundance per site and method
Transect_mean <- Transect_long %>%
      group_by(Site, Latitude, Species) %>%
     summarise(Mean_abundance = mean(Abundance),n = n(), sd = sd(Abundance), se = sd/sqrt(n))
  • Calculate abundances by group

Follow the same steps as abundances by type

TransectGroup_sum <- Transect_metadata %>%
  group_by(Site, Latitude, Census, Group) %>%
  summarise(Total_abundance = sum(Individuals))

TransectGroup_wide <- TransectGroup_sum %>%
  ungroup()  %>%
  pivot_wider(names_from = Group, values_from = Total_abundance,
              values_fill = 0) 
TransectGroup_wide <- TransectGroup_wide%>%
  mutate(Herbivores = Herbivore_fish + Herbivore_gastropod + Urchin + Omnivorous_crab, Invertebrate_predators = Predatory_gastropod+ Snail + Predatory_crab + Sea_star + Omnivorous_crab + Schrimp)


TransectGroup_long <- TransectGroup_wide %>%
  gather(Group, Abundance, Herbivore_fish:Invertebrate_predators)

TransectGroup_mean <- TransectGroup_long %>%
  group_by(Site, Latitude, Group) %>%
  summarise(Total_abundance = mean(Abundance),n = n(), sd = sd(Abundance), se = sd/sqrt(n))
  1. Calculate Species Richness

First, add a column of presence when abundance is higher than 0.

Transect_long <- Transect_long %>%
  mutate(Presence = if_else(Abundance > 0,1,0))

Calculate sum of all species in the same site and census.

Richness_transect <- Transect_long %>%
  group_by(Site, Latitude, Census) %>%
  summarise(Richness = sum(Presence))

#Calculate mean by site  
Richness_mean <- Richness_transect %>%
  group_by(Site, Latitude)%>%
  summarise(Mean_SR = mean(Richness),n = n(), sd = sd(Richness), se = sd/sqrt(n))
  1. Add abundances to consumption data

Calculate means and change it to the wide format

#Per Group
Group_abundances_mean <- TransectGroup_long%>%
  filter(Group == "Invertebrate_predators" | Group == "Predatory_fish"| Group == "Herbivores")%>%
  group_by(Site, Group)%>%
  summarise(Mean_Abundance = mean(Abundance)) %>%
  pivot_wider(names_from = Group, values_from = Mean_Abundance,
              values_fill = 0)

Then, add abundances data to consumption

  • Predation_env + group abundance summarised data
#Select only the needed columns
Group_abundances_mean_pred <- dplyr::select(Group_abundances_mean, Site, Invertebrate_predators, Predatory_fish)

Predation_4 <- Predation_3 %>%
    left_join(Group_abundances_mean_pred, by = "Site")
  • Herbivory_env + Benthic summarised data
#Select only the needed columns
Group_abundances_mean_herb <- dplyr::select(Group_abundances_mean, Site, Herbivores)

Herbivory_3 <- Herbivory_2 %>%
  left_join(Group_abundances_mean_herb, by = "Site")
  1. Add species richness to consumption data

Calculate means and change it to the wide format

Richness_meanpred <- Richness_transect %>%
  group_by(Site)%>%
  summarise(Mean_SR = mean(Richness)) 

Then, add abundances data to consumption

  • Predation_4 + species richness summarised data
Predation_5 <- Predation_4 %>%
    left_join(Richness_meanpred, by = "Site")
  • Herbivory_3 + species richness summarised data
Herbivory_4 <- Herbivory_3 %>%
  left_join(Richness_meanpred, by = "Site")

2e): Community composition

Summarise abundances of species per site and standardise

Transect_sum <- Transect_long %>%
  group_by(Site, Latitude, Species) %>%
  summarise(Total_abundance = mean(Abundance))
## `summarise()` has grouped output by 'Site', 'Latitude'. You can override using
## the `.groups` argument.
Transect_wide2 <- Transect_sum %>%
ungroup() %>%
pivot_wider(names_from = Species, values_from = Total_abundance, values_fill = 0)

Transect_val <- Transect_wide2 %>%
  dplyr::select(- Site, -Latitude) 

Transect_stand <- decostand(Transect_val, method = "total")

Calculate Bray-Curtis dissimilariry

Transect_comp <- vegdist(Transect_stand, "bray")

Run Principal Coordinate Analysis (PcoA) to obtain the component explaining the highers variation

pcoa <- cmdscale(Transect_comp, eig = TRUE, add = TRUE)
pcoa$eig
##  [1]  1.740530e+00  1.329554e+00  8.438524e-01  6.203963e-01  5.740214e-01
##  [6]  3.517972e-01  1.671264e-01  7.777519e-02  5.878784e-02  3.878297e-02
## [11] -7.537429e-17 -2.615278e-16
(pcoa$eig/sum(pcoa$eig))*100 #proportion of variation
##  [1]  2.999557e+01  2.291298e+01  1.454260e+01  1.069165e+01  9.892446e+00
##  [6]  6.062727e+00  2.880187e+00  1.340345e+00  1.013125e+00  6.683695e-01
## [11] -1.298969e-15 -4.507061e-15
#convert pcoa results into data frame that can be plotted
pcoa_df <- data.frame(pcoa$points)
colnames(pcoa_df) <- c("PCo1", "PCo2")
pcoa_df <- pcoa_df %>%
  dplyr::select(-c(PCo2))
pcoa_df$Site <- factor(Transect_wide2$Site)
pcoa_df$Latitude <- Transect_wide2$Latitude

PCo 1 represents 30% of the variation community composition

  1. Add PCos to consumption data
  • Predation_4 + PCo
pcoa_df1<- dplyr::select(pcoa_df, -Latitude)
Predation_6 <- Predation_5 %>%
    left_join(pcoa_df1, by = "Site")
  • Herbivory_3 + PCo
Herbivory_5 <- Herbivory_4  %>%
    left_join(pcoa_df1, by = "Site")

2f): Phlorotannins

  1. Add metadata
Phlorotannins_metadata <- Phlorotannins %>%
left_join(Metadata, by = c("Site"))
  1. Summarise phlorotannins per site
Phlorotannins_mean <- Phlorotannins_metadata %>%
  group_by(Site, Latitude) %>%
  summarise(phloro_mean = mean(Phlorotannins), n = n(), sd = sd(Phlorotannins), se = sd/sqrt(n))
  1. Herbivory + phlorotannins
Phlorotannins_sum <- Phlorotannins %>%
  group_by(Site) %>%
  summarise(phloro_mean = mean(Phlorotannins))
Herbivory_6 <- Herbivory_5  %>%
    left_join(Phlorotannins_sum, by = "Site")

Part 3: Latitudinal patterns in explanatory variables

Latitudinal patterns were analysed for the recorded explanatory variables using linear models.

3a): Explore latitudinal relationships

Set general theme for plots

Theme1 <- theme_bw() + 
    theme(panel.grid.major = element_line(colour = "white"),
          panel.grid.minor = element_line(colour = "white"), axis.text.y = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 9), axis.title.x =     element_text(size = 9))

For these analyses we used the approach and code of model selection for environmental and spatial gradients provided by Anderson et al. (2022)

1) Temperature vs latitude

a. Analyses

For future plotting

Environmental_2$Latitude <- Environmental_2$Latitude*-1
#change sign of latitude 
# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$Temperature, xlab = "Latitude", ylab = "Temperature", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$Temperature, pch = 21, col = mygrey, bg = mygrey)

descdist(Environmental_2$Temperature, discrete = FALSE)

## summary statistics
## ------
## min:  10.2   max:  18.7 
## median:  16 
## mean:  15.32917 
## estimated sd:  2.365349 
## estimated skewness:  -0.7647711 
## estimated kurtosis:  2.651633
normal_dist <- fitdist(Environmental_2$Temperature, "weibull")
plot(normal_dist)

T_mod1 <- glmer(Temperature ~ Latitude+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))
T_mod2 <- glmer(Temperature ~poly(Latitude,2)+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))
T_mod3 <- glmer(Temperature ~poly(Latitude,3)+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))
AICc(T_mod1, T_mod3, T_mod2)
##        df     AICc
## T_mod1  4 80.68285
## T_mod3  6 84.78450
## T_mod2  5 82.55903
r2(T_mod1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.930
##      Marginal R2: 0.713
r2(T_mod2)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.935
##      Marginal R2: 0.744
r2(T_mod3)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.937
##      Marginal R2: 0.761

T_mod1 is the best fitted model

res <- resid(T_mod1)
#produce residual vs. fitted plot
plot(fitted(T_mod1), res)

#add a horizontal line at 0 
abline(0,0)

#create Q-Q plot for residuals
qqnorm(res)

#add a straight diagonal line to the plot
qqline(res)

summary(T_mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: Temperature ~ Latitude + (1 | Site)
##    Data: Environmental_2
## 
##      AIC      BIC   logLik deviance df.resid 
##     79.8     87.2    -35.9     71.8       44 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.34203 -0.41973  0.05954  0.42096  1.54302 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 0.004483 0.06696 
##  Residual             0.001452 0.03811 
## Number of obs: 48, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error t value Pr(>|z|)    
## (Intercept)  3.23945    0.20465  15.829  < 2e-16 ***
## Latitude    -0.01755    0.00669  -2.623  0.00871 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## Latitude -0.975
r2(T_mod1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.930
##      Marginal R2: 0.713

b. Plot predicted relationship

predictions <- ggpredict(T_mod1, terms ="Latitude")
# Plot the GAM model predictions
Temperature_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=Temperature), width = 0.08, height = 0.03, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.high, ymax = conf.low), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y =predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Temperature (°C)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

2) Salinity vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$Salinity_ppt, xlab = "Latitude", ylab = "Salinity", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$Salinity_ppt, pch = 21, col = mygrey, bg = mygrey)

descdist(Environmental_2$Salinity_ppt, discrete = FALSE)

## summary statistics
## ------
## min:  32.8   max:  34.4 
## median:  33.8 
## mean:  33.73125 
## estimated sd:  0.3882071 
## estimated skewness:  -0.9322086 
## estimated kurtosis:  3.328625
normal_dist <- fitdist(Environmental_2$Salinity_ppt, "weibull")
plot(normal_dist)

S_mod1 <- glmer(Salinity_ppt ~ Latitude+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0215438 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
S_mod2 <- glmer(Salinity_ppt ~poly(Latitude,2)+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))
S_mod3 <- glmer(Salinity_ppt ~poly(Latitude,3)+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))
check_model(S_mod1)

AICc(S_mod1, S_mod3, S_mod2)
##        df     AICc
## S_mod1  4 20.65664
## S_mod3  6 24.12187
## S_mod2  5 21.50425
BIC(S_mod1, S_mod3, S_mod2)
##        df      BIC
## S_mod1  4 27.21121
## S_mod3  6 33.30029
## S_mod2  5 29.43168
summary(S_mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: Salinity_ppt ~ Latitude + (1 | Site)
##    Data: Environmental_2
## 
##      AIC      BIC   logLik deviance df.resid 
##     19.7     27.2     -5.9     11.7       44 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8003 -0.5564  0.1690  0.5914  1.4235 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Site     (Intercept) 2.189e-05 0.004679
##  Residual             5.926e-05 0.007698
## Number of obs: 48, groups:  Site, 12
## 
## Fixed effects:
##               Estimate Std. Error t value Pr(>|z|)    
## (Intercept)  3.5487577  0.0083718 423.896  < 2e-16 ***
## Latitude    -0.0010191  0.0002757  -3.696 0.000219 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## Latitude -0.953
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0215438 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
r2(S_mod1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.547
##      Marginal R2: 0.379

b. Plot predicted relationship

predictions <- ggpredict(S_mod1, terms ="Latitude [all]")
# Plot the GAM model predictions
Salinity_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=Salinity_ppt), width = 0.08, height = 0.05, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data =predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Salinity (ppt)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

3) Oxygen vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$D_oxygen_mgL, xlab = "Latitude", ylab = "DO", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$D_oxygen_mgL, pch = 21, col = mygrey, bg = mygrey)

descdist(Environmental_2$D_oxygen_mgL, discrete = FALSE)

## summary statistics
## ------
## min:  1.62   max:  12.15 
## median:  6.445 
## mean:  6.381042 
## estimated sd:  2.626742 
## estimated skewness:  0.3329725 
## estimated kurtosis:  2.893175
normal_dist <- fitdist(Environmental_2$D_oxygen_mgL, "norm")
plot(normal_dist)

Ox_mod1 <- glmer(D_oxygen_mgL ~ Latitude+ (1|Site),  data= Environmental_2, family = gaussian)
## Warning in glmer(D_oxygen_mgL ~ Latitude + (1 | Site), data = Environmental_2,
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly
Ox_mod2 <- glmer(D_oxygen_mgL ~poly(Latitude,2)+ (1|Site),  data= Environmental_2, family = gaussian)
## Warning in glmer(D_oxygen_mgL ~ poly(Latitude, 2) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
Ox_mod3 <- glmer(D_oxygen_mgL ~poly(Latitude,3)+ (1|Site),  data= Environmental_2, family = gaussian)
## Warning in glmer(D_oxygen_mgL ~ poly(Latitude, 3) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
AICc(Ox_mod1, Ox_mod2, Ox_mod3)
##         df     AICc
## Ox_mod1  4 221.6313
## Ox_mod2  5 209.4476
## Ox_mod3  6 207.2213
BIC(Ox_mod1, Ox_mod2, Ox_mod3)
##         df      BIC
## Ox_mod1  4 228.1858
## Ox_mod2  5 217.3750
## Ox_mod3  6 216.3998

Improvement cubic model shows the best fit

Two options may be possible based on AICc. Very low df in Ox_mod 3

summary(Ox_mod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: D_oxygen_mgL ~ poly(Latitude, 3) + (1 | Site)
##    Data: Environmental_2
## 
## REML criterion at convergence: 193.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.37338 -0.33408  0.05551  0.26713  3.04460 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 4.611    2.147   
##  Residual             3.029    1.741   
## Number of obs: 48, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)          6.3810     0.6688   9.540
## poly(Latitude, 3)1   0.6637     4.6339   0.143
## poly(Latitude, 3)2  -6.5554     4.6339  -1.415
## poly(Latitude, 3)3  -0.1882     4.6339  -0.041
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 0.000                 
## ply(Ltt,3)2 0.000  0.000          
## ply(Ltt,3)3 0.000  0.000   0.000
#Plot the GAM model predictions
Oxygen_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=D_oxygen_mgL), width = 0.08, height = 0.07, size = 1, alpha = 0.6)  + 
  Theme1 +
labs( y="Dissolved oxygen (mg"~L^-1*")", x="Latitude (°S)") +
theme(axis.title.x=element_text(size = 10),
        axis.text.x=element_text(size = 10), axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 8))
par(mfrow = c(2, 2))
check_model(Ox_mod3)

4) Visibility vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$V_Visibility, xlab = "Latitude", ylab = "Visibility", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$V_Visibility, pch = 21, col = mygrey, bg = mygrey)

descdist(Environmental_2$V_Visibility, discrete = FALSE)

## summary statistics
## ------
## min:  3   max:  10 
## median:  5.75 
## mean:  5.854167 
## estimated sd:  1.994563 
## estimated skewness:  0.4867126 
## estimated kurtosis:  2.512658
normal_dist <- fitdist(Environmental_2$V_Visibility, "norm")
plot(normal_dist)

V_mod1 <- glmer(V_Visibility ~ Latitude+ (1|Site),  data= Environmental_2, family = gaussian())
## Warning in glmer(V_Visibility ~ Latitude + (1 | Site), data = Environmental_2,
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly
V_mod2 <- glmer(V_Visibility ~poly(Latitude,2)+ (1|Site),  data= Environmental_2, family = gaussian())
## Warning in glmer(V_Visibility ~ poly(Latitude, 2) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
V_mod3 <- glmer(V_Visibility ~poly(Latitude,3)+ (1|Site),  data= Environmental_2, family = gaussian())
## Warning in glmer(V_Visibility ~ poly(Latitude, 3) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
check_model(V_mod1)

AICc(V_mod1, V_mod2, V_mod3)
##        df     AICc
## V_mod1  4 171.7591
## V_mod2  5 162.0855
## V_mod3  6 158.9580
BIC(V_mod1, V_mod2, V_mod3)
##        df      BIC
## V_mod1  4 178.3137
## V_mod2  5 170.0129
## V_mod3  6 168.1365

V_mod3 is the best fitted

summary(V_mod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: V_Visibility ~ poly(Latitude, 3) + (1 | Site)
##    Data: Environmental_2
## 
## REML criterion at convergence: 144.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.83174 -0.54808  0.07988  0.44641  2.45295 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 3.2862   1.8128  
##  Residual             0.8715   0.9336  
## Number of obs: 48, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)          5.8542     0.5404  10.834
## poly(Latitude, 3)1   5.0601     3.7438   1.352
## poly(Latitude, 3)2  -0.1325     3.7438  -0.035
## poly(Latitude, 3)3  -4.2251     3.7438  -1.129
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 0.000                 
## ply(Ltt,3)2 0.000  0.000          
## ply(Ltt,3)3 0.000  0.000   0.000
# Plot the GAM model predictions
Visibility_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=V_Visibility), width = 0.08, height = 0.08, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="Visibility (m)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

5) PCOA vs latitude

a. Analyses

pcoa_df$Latitude <- pcoa_df$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
   plot(pcoa_df$Latitude, pcoa_df$PCo1, xlab = "Latitude", ylab = "PCoA1", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(pcoa_df$Latitude, pcoa_df$PCo1, pch = 21, col = mygrey, bg = mygrey)

descdist(pcoa_df$PCo1, discrete = FALSE)

## summary statistics
## ------
## min:  -0.3828236   max:  0.6887275 
## median:  -0.1217363 
## mean:  5.247539e-17 
## estimated sd:  0.3977814 
## estimated skewness:  0.9122292 
## estimated kurtosis:  2.331104
normal_dist <- fitdist(log(pcoa_df$PCo1+1), "norm")
plot(normal_dist)

PCO_mod1 <- glm(log(PCo1 +1) ~ Latitude,  data= pcoa_df, family = gaussian())
PCO_mod2 <- glm(log(PCo1 +1) ~poly(Latitude,2),  data= pcoa_df, family = gaussian())
PCO_mod3 <- glm(log(PCo1 +1) ~poly(Latitude,3),  data= pcoa_df, family = gaussian())
check_model(PCO_mod2)

check_model(PCO_mod1)

AICc(PCO_mod1,PCO_mod2, PCO_mod3)
##          df     AICc
## PCO_mod1  3 14.72754
## PCO_mod2  4 17.28033
## PCO_mod3  5 23.21428
summary(PCO_mod1)
## 
## Call:
## glm(formula = log(PCo1 + 1) ~ Latitude, family = gaussian(), 
##     data = pcoa_df)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.85098    0.43453  -1.958   0.0786 .
## Latitude     0.02632    0.01421   1.853   0.0936 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1132384)
## 
##     Null deviance: 1.5211  on 11  degrees of freedom
## Residual deviance: 1.1324  on 10  degrees of freedom
## AIC: 11.728
## 
## Number of Fisher Scoring iterations: 2
par(mfrow = c(2, 2))
gam.check(PCO_mod1)

## 
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.

No significant effect

# Plot the GAM model predictions
PCOA_plot <- ggplot() +
  geom_jitter(data = pcoa_df, aes(x= Latitude, y=PCo1), width = 0, height = 0, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="log(PCoA1 + 1)", x="Latitude (°S)")

6) Kelp vs latitude

a. Analyses

For further plotting

Benthic_metadata_wide$Latitude <- Benthic_metadata_wide$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
   plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Kelp, xlab = "Latitude", ylab = "Kelp", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Kelp, pch = 21, col = mygrey, bg = mygrey)

Kelp_mod1 <- glmmTMB(Kelp ~ Latitude+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
Kelp_mod2 <- glmmTMB(Kelp ~poly(Latitude,2)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
Kelp_mod3 <- glmmTMB(Kelp ~poly(Latitude,3)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
check_model(Kelp_mod3)
## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(Kelp_mod2, Kelp_mod1, Kelp_mod3)
##           df     AICc
## Kelp_mod2  6 470.8080
## Kelp_mod1  5 468.5267
## Kelp_mod3  7 473.2009
BIC(Kelp_mod2, Kelp_mod1, Kelp_mod3)
##           df      BIC
## Kelp_mod2  6 483.7652
## Kelp_mod1  5 479.4769
## Kelp_mod3  7 488.0979
summary(Kelp_mod1)
##  Family: tweedie  ( log )
## Formula:          Kelp ~ Latitude + (1 | Site)
## Data: Benthic_metadata_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    467.7    479.5   -228.8    457.7       73 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 1.159    1.076   
## Number of obs: 78, groups:  Site, 12
## 
## Dispersion parameter for tweedie family (): 6.27 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  6.42602    1.70405   3.771 0.000163 ***
## Latitude    -0.15074    0.06016  -2.506 0.012217 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(Kelp_mod1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.800
##      Marginal R2: 0.371

b. Plot predicted relationship

predictions <- ggpredict(Kelp_mod1, terms ="Latitude [all]")
# Plot the GAM model predictions
Kelp_plot <- ggplot() +
  geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Kelp), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Kelp cover (%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y=element_text(8))

7) Bare_rock vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Bare_rock, xlab = "Latitude", ylab = "Bare_rock", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Bare_rock, pch = 21, col = mygrey, bg = mygrey)

BR_mod1 <- mgcv::gam(Bare_rock ~ s(Latitude), gamma= 1.4, data = Benthic_metadata_wide, family = gaussian())
BR_mod2 <- mgcv::gam(Bare_rock ~ s(Latitude), gamma=1.4, data = Benthic_metadata_wide, family = tw())


descdist(Benthic_metadata_wide$Bare_rock, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  100 
## median:  23.85013 
## mean:  36.00036 
## estimated sd:  32.50373 
## estimated skewness:  0.5844581 
## estimated kurtosis:  1.884422
normal_dist <- fitdist(Benthic_metadata_wide$Bare_rock, "exp")
plot(normal_dist)

BR_mod1 <- glmmTMB(Bare_rock ~ Latitude+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
BR_mod2 <- glmmTMB(Bare_rock ~poly(Latitude,2)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
BR_mod3 <- glmmTMB(Bare_rock ~poly(Latitude,3)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
check_model(BR_mod2)
## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(BR_mod1, BR_mod2, BR_mod3)
##         df     AICc
## BR_mod1  5 634.5401
## BR_mod2  6 634.5044
## BR_mod3  7 636.9098
BIC(BR_mod1, BR_mod2, BR_mod3)
##         df      BIC
## BR_mod1  5 645.4903
## BR_mod2  6 647.4616
## BR_mod3  7 651.8068
summary(BR_mod2)
##  Family: tweedie  ( log )
## Formula:          Bare_rock ~ poly(Latitude, 2) + (1 | Site)
## Data: Benthic_metadata_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    633.3    647.5   -310.7    621.3       72 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.8647   0.9299  
## Number of obs: 78, groups:  Site, 12
## 
## Dispersion parameter for tweedie family ():  4.1 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           3.126      0.283  11.046   <2e-16 ***
## poly(Latitude, 2)1    2.890      2.449   1.180    0.238    
## poly(Latitude, 2)2   -3.909      2.433  -1.606    0.108    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(BR_mod2)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.859
##      Marginal R2: 0.225

R2 of mod2 is higher

# Plot the GAM model predictions
Bare_plot <- ggplot() +
  geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Bare_rock), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="Bare rock cover(%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y=element_text(8))

8) Understorey_algae vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Understorey_algae, xlab = "Latitude", ylab = "Understorey_algae", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Understorey_algae, pch = 21, col = mygrey, bg = mygrey)

descdist(Benthic_metadata_wide$Understorey_algae, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  82.92683 
## median:  0 
## mean:  10.92887 
## estimated sd:  19.35412 
## estimated skewness:  2.110129 
## estimated kurtosis:  7.04205
normal_dist <- fitdist(Benthic_metadata_wide$Understorey_algae, "exp")
plot(normal_dist)

UA_mod1 <- glmmTMB(Understorey_algae ~ Latitude+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
UA_mod2 <- glmmTMB(Understorey_algae ~poly(Latitude,2)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
UA_mod3 <- glmmTMB(Understorey_algae ~poly(Latitude,3)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
check_model(UA_mod2)
## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(UA_mod1, UA_mod2, UA_mod3)
##         df     AICc
## UA_mod1  5 342.5629
## UA_mod2  6 344.5077
## UA_mod3  7 342.1709
BIC(UA_mod1, UA_mod2, UA_mod3)
##         df      BIC
## UA_mod1  5 353.5131
## UA_mod2  6 357.4649
## UA_mod3  7 357.0679
summary(UA_mod1)
##  Family: tweedie  ( log )
## Formula:          Understorey_algae ~ Latitude + (1 | Site)
## Data: Benthic_metadata_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    341.7    353.5   -165.9    331.7       73 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 5.819    2.412   
## Number of obs: 78, groups:  Site, 12
## 
## Dispersion parameter for tweedie family (): 6.61 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  1.75831    3.33416   0.527    0.598
## Latitude    -0.04022    0.11085  -0.363    0.717
r2(UA_mod1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.934
##      Marginal R2: 0.011
# Plot the GAM model predictions
UA_plot <- ggplot() +
  geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Understorey_algae), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="UA cover(%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y=element_text(8))

9) Predatory fish abundance vs latitude

TransectGroup_wide$Latitude <- TransectGroup_wide$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
   plot(TransectGroup_wide$Latitude, TransectGroup_wide$Predatory_fish, xlab = "Latitude", ylab = "PF", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(TransectGroup_wide$Latitude, TransectGroup_wide$Predatory_fish, pch = 21, col = mygrey, bg = mygrey)

descdist(TransectGroup_wide$Predatory_fish, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  333 
## median:  29 
## mean:  51.61538 
## estimated sd:  75.31432 
## estimated skewness:  2.720417 
## estimated kurtosis:  10.84418
normal_dist <- fitdist(TransectGroup_wide$Predatory_fish, "exp")
plot(normal_dist)

PF_mod1 <- glmer.nb(Predatory_fish ~ Latitude+ (1|Site),  data= TransectGroup_wide)
PF_mod2 <- glmer.nb(Predatory_fish ~poly(Latitude,2)+ (1|Site),  data= TransectGroup_wide)
PF_mod3 <- glmer.nb(Predatory_fish ~poly(Latitude,3)+ (1|Site),  data= TransectGroup_wide)
## boundary (singular) fit: see help('isSingular')
check_model(PF_mod3)

AICc(PF_mod2, PF_mod1, PF_mod3)
##         df     AICc
## PF_mod2  5 257.1651
## PF_mod1  4 257.3351
## PF_mod3  6 249.9318
BIC(PF_mod2, PF_mod1, PF_mod3)
##         df      BIC
## PF_mod2  5 260.4556
## PF_mod1  4 260.4627
## PF_mod3  6 253.0594

mod 3 is singular, so mod 1 was chosen

summary(PF_mod2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(1.2985)  ( log )
## Formula: Predatory_fish ~ poly(Latitude, 2) + (1 | Site)
##    Data: TransectGroup_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    254.2    260.5   -122.1    244.2       21 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0949 -0.6533 -0.2637  0.1770  2.9381 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4711   0.6864  
## Number of obs: 26, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          3.3709     0.3024  11.147  < 2e-16 ***
## poly(Latitude, 2)1  -4.8871     1.5132  -3.230  0.00124 ** 
## poly(Latitude, 2)2  -2.7915     1.5243  -1.831  0.06705 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,2)1
## ply(Ltt,2)1 0.071         
## ply(Ltt,2)2 0.159  0.175
r2(PF_mod2)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.749
##      Marginal R2: 0.546

b. Plot predicted relationship

predictions <- ggpredict(PF_mod2, terms ="Latitude [all]")
# Plot the GAM model predictions
PredFish_plot <- ggplot() +
  geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Predatory_fish), width = 0.12, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Predatory Fish", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y=element_text(8))

10) Herbivorous fish abundance vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(TransectGroup_wide$Latitude, TransectGroup_wide$Herbivores, xlab = "Latitude", ylab = "Herbivores", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(TransectGroup_wide$Latitude, TransectGroup_wide$Herbivores, pch = 21, col = mygrey, bg = mygrey)

descdist(TransectGroup_wide$Herbivores, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  1912 
## median:  155.5 
## mean:  389.4231 
## estimated sd:  547.3004 
## estimated skewness:  1.702098 
## estimated kurtosis:  5.106245
normal_dist <- fitdist(TransectGroup_wide$Herbivores, "exp")
plot(normal_dist)

Herb_mod1 <- glmer.nb(Herbivores ~ Latitude+ (1|Site),  data= TransectGroup_wide)
Herb_mod2 <- glmer.nb(Herbivores ~poly(Latitude,2)+ (1|Site),  data= TransectGroup_wide)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0495445 (tol = 0.002, component 1)
Herb_mod3 <- glmer.nb(Herbivores ~poly(Latitude,3)+ (1|Site),  data= TransectGroup_wide)
check_model(Herb_mod1)

AICc(Herb_mod1, Herb_mod2, Herb_mod3)
##           df     AICc
## Herb_mod1  4 326.4081
## Herb_mod2  5 326.9362
## Herb_mod3  6 327.9244
BIC(Herb_mod1, Herb_mod2, Herb_mod3)
##           df      BIC
## Herb_mod1  4 329.5358
## Herb_mod2  5 330.2266
## Herb_mod3  6 331.0519
summary(Herb_mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(1.5227)  ( log )
## Formula: Herbivores ~ Latitude + (1 | Site)
##    Data: TransectGroup_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    324.5    329.5   -158.3    316.5       22 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2022 -0.4869 -0.0675  0.3265  1.3791 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 5.748    2.398   
## Number of obs: 26, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  10.1339     3.3117   3.060  0.00221 **
## Latitude     -0.2029     0.1101  -1.844  0.06523 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## Latitude -0.976
r2(Herb_mod1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.939
##      Marginal R2: 0.250
# Plot the GAM model predictions
Herb_plot <- ggplot() +
  geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Herbivores), width = 0.2, height = 0.5, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="Herbivores", x="Latitude (°S)") 

11) Invertebrate predators abundance vs latitude

# First, examine a simple scatter plot (always a good idea!)
   plot(TransectGroup_wide$Latitude, TransectGroup_wide$Invertebrate_predators, xlab = "Latitude", ylab = "Invertebrate pred", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(TransectGroup_wide$Latitude, TransectGroup_wide$Invertebrate_predators, pch = 21, col = mygrey, bg = mygrey)

descdist(TransectGroup_wide$Invertebrate_predators, discrete = FALSE)

## summary statistics
## ------
## min:  8   max:  753 
## median:  174.5 
## mean:  242.3077 
## estimated sd:  208.9332 
## estimated skewness:  1.010871 
## estimated kurtosis:  3.417233
normal_dist <- fitdist(TransectGroup_wide$Invertebrate_predators, "exp")
plot(normal_dist)

Inv_mod1 <- glmer.nb(Invertebrate_predators ~ Latitude+ (1|Site),  data= TransectGroup_wide)
Inv_mod2 <- glmer.nb(Invertebrate_predators ~poly(Latitude,2)+ (1|Site),  data= TransectGroup_wide)
Inv_mod3 <- glmer.nb(Invertebrate_predators ~poly(Latitude,3)+ (1|Site),  data= TransectGroup_wide)
check_model(Inv_mod1)

AICc(Inv_mod2, Inv_mod1, Inv_mod3)
##          df     AICc
## Inv_mod2  5 337.9285
## Inv_mod1  4 337.7932
## Inv_mod3  6 337.1040
BIC(Inv_mod2, Inv_mod1, Inv_mod3)
##          df      BIC
## Inv_mod2  5 341.2190
## Inv_mod1  4 340.9208
## Inv_mod3  6 340.2315
summary(Inv_mod3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(3.6522)  ( log )
## Formula: Invertebrate_predators ~ poly(Latitude, 3) + (1 | Site)
##    Data: TransectGroup_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    332.7    340.2   -160.3    320.7       20 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.48939 -0.60016 -0.01322  0.64400  1.52814 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4727   0.6875  
## Number of obs: 26, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          5.0982     0.2268  22.483   <2e-16 ***
## poly(Latitude, 3)1  -1.1682     1.1803  -0.990   0.3223    
## poly(Latitude, 3)2  -2.3212     1.1633  -1.995   0.0460 *  
## poly(Latitude, 3)3  -2.6746     1.1575  -2.311   0.0208 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 -0.067                
## ply(Ltt,3)2  0.072 -0.069         
## ply(Ltt,3)3 -0.017  0.024  -0.085
r2(Inv_mod3)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.808
##      Marginal R2: 0.437
predictions <- ggpredict(Inv_mod3, terms ="Latitude [all]")
# Plot the GAM model predictions
Inv_plot <- ggplot() +
  geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Invertebrate_predators), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
    geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Invertebrate predators", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y = element_text(size = 7))

12) Species richness vs latitude

Richness_transect$Latitude <- Richness_transect$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
   plot(Richness_transect$Latitude, Richness_transect$Richness, xlab = "Latitude", ylab = "Richness", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Richness_transect$Latitude, Richness_transect$Richness, pch = 21, col = mygrey, bg = mygrey)

descdist(Richness_transect$Richness, discrete = FALSE)

## summary statistics
## ------
## min:  2   max:  15 
## median:  8 
## mean:  8.115385 
## estimated sd:  3.076712 
## estimated skewness:  0.3378056 
## estimated kurtosis:  3.339398
normal_dist <- fitdist(Richness_transect$Richness, "norm")
plot(normal_dist)

R_mod1 <- glmer(Richness ~ Latitude+ (1|Site),  data= Richness_transect, family = gaussian())
## Warning in glmer(Richness ~ Latitude + (1 | Site), data = Richness_transect, :
## calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly
R_mod2 <- glmer(Richness ~poly(Latitude,2)+ (1|Site),  data= Richness_transect, family = gaussian())
## Warning in glmer(Richness ~ poly(Latitude, 2) + (1 | Site), data =
## Richness_transect, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
R_mod3 <- glmer(Richness ~poly(Latitude,3)+ (1|Site),  data= Richness_transect, family = gaussian())
## Warning in glmer(Richness ~ poly(Latitude, 3) + (1 | Site), data =
## Richness_transect, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
AICc(R_mod2, R_mod1, R_mod3)
##        df     AICc
## R_mod2  5 127.3899
## R_mod1  4 136.8463
## R_mod3  6 124.1031
BIC(R_mod2, R_mod1, R_mod3)
##        df      BIC
## R_mod2  5 130.6804
## R_mod1  4 139.9739
## R_mod3  6 127.2306
summary(R_mod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Richness ~ poly(Latitude, 3) + (1 | Site)
##    Data: Richness_transect
## 
## REML criterion at convergence: 107.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6252 -0.5226  0.0082  0.3615  2.0845 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 5.464    2.338   
##  Residual             4.171    2.042   
## Number of obs: 26, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)          8.0744     0.7901  10.220
## poly(Latitude, 3)1  -3.4362     4.0846  -0.841
## poly(Latitude, 3)2  -3.0754     4.0715  -0.755
## poly(Latitude, 3)3  -5.9172     4.0749  -1.452
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 -0.079                
## ply(Ltt,3)2  0.057 -0.081         
## ply(Ltt,3)3 -0.014  0.039  -0.072
r2(R_mod3)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.649
##      Marginal R2: 0.189

No significant effect

# Plot the GAM model predictions
Richness_plot <- ggplot() +
  geom_jitter(data = Richness_transect, aes(x= Latitude, y=Richness), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="Species Richness", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y = element_text(size = 8))

13) Fish length vs latitude

Size_fish <- filter(Size_mean_sp, Type == "Fish")
Size_fish$Latitude <- Size_fish$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
   plot(Size_fish$Latitude, Size_fish$mean_Size_sp, xlab = "Latitude", ylab = "Fish Length", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Size_fish$Latitude, Size_fish$mean_Size_sp, pch = 21, col = mygrey, bg = mygrey)

descdist(Size_fish$mean_Size_sp, discrete = FALSE)

## summary statistics
## ------
## min:  5.4   max:  65 
## median:  23 
## mean:  22.49199 
## estimated sd:  9.684538 
## estimated skewness:  0.9708081 
## estimated kurtosis:  5.717812
normal_dist <- fitdist(Size_fish$mean_Size_sp, "logis")
plot(normal_dist)

Length_mod1 <- glmer(mean_Size_sp ~ Latitude+ (1|Site) + (1|Species),  data= Size_fish, family = Gamma(link = "log"))
Length_mod2 <- glmer(mean_Size_sp ~poly(Latitude,2)+ (1|Site) + (1|Species),  data= Size_fish, family = Gamma(link = "log"))
Length_mod3 <- glmer(mean_Size_sp ~poly(Latitude,3)+ (1|Site) + (1|Species),  data= Size_fish, family = Gamma(link = "log"))
check_model(Length_mod1)

AICc(Length_mod1, Length_mod2, Length_mod3)
##             df     AICc
## Length_mod1  5 773.3847
## Length_mod2  6 772.3818
## Length_mod3  7 774.6382
BIC(Length_mod1, Length_mod2, Length_mod3)
##             df      BIC
## Length_mod1  5 786.7493
## Length_mod2  6 788.3065
## Length_mod3  7 793.0830
summary(Length_mod2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: mean_Size_sp ~ poly(Latitude, 2) + (1 | Site) + (1 | Species)
##    Data: Size_fish
## 
##      AIC      BIC   logLik deviance df.resid 
##    771.6    788.3   -379.8    759.6      113 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2255 -0.4292  0.0110  0.4239  4.0786 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Species  (Intercept) 0.091860 0.30308 
##  Site     (Intercept) 0.002929 0.05412 
##  Residual             0.073020 0.27022 
## Number of obs: 119, groups:  Species, 28; Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value Pr(>|z|)    
## (Intercept)          3.0014     0.1174  25.566   <2e-16 ***
## poly(Latitude, 2)1   0.1303     0.3940   0.331   0.7409    
## poly(Latitude, 2)2  -0.7157     0.3493  -2.049   0.0405 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,2)1
## ply(Ltt,2)1 -0.085        
## ply(Ltt,2)2 -0.027  0.036
r2(Length_mod2)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.585
##      Marginal R2: 0.026
predictions <- ggpredict(Length_mod2, terms ="Latitude [all]")
# Plot the GAM model predictions
Size_plot <- ggplot() +
  geom_jitter(data = Size_fish, aes(x= Latitude, y=mean_Size_sp), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Fish length (cm)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y=element_text(8))

14) Phlorotannins vs latitude

Phlorotannins_metadata$Latitude <- Phlorotannins_metadata$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
   plot(Phlorotannins_metadata$Latitude, Phlorotannins_metadata$Phlorotannins, xlab = "Latitude", ylab = "phlorotannins", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Phlorotannins_metadata$Latitude, Phlorotannins_metadata$Phlorotannins, pch = 21, col = mygrey, bg = mygrey)

descdist(Phlorotannins_metadata$Phlorotannins, discrete = FALSE)

## summary statistics
## ------
## min:  1.83   max:  94.69 
## median:  11.625 
## mean:  17.26067 
## estimated sd:  15.2777 
## estimated skewness:  2.275762 
## estimated kurtosis:  9.890266
normal_dist <- fitdist(Phlorotannins_metadata$Phlorotannins, "exp")
plot(normal_dist)

Phloro_mod1 <- glmer(Phlorotannins ~ Latitude+ (1|Site),  data= Phlorotannins_metadata, family = Gamma(link= "log"))
Phloro_mod2 <- glmer(Phlorotannins ~poly(Latitude,2)+ (1|Site),  data= Phlorotannins_metadata, family = Gamma(link= "log"))
Phloro_mod3 <- glmer(Phlorotannins ~poly(Latitude,3)+ (1|Site),  data= Phlorotannins_metadata, family = Gamma(link= "log"))
check_model(Phloro_mod1)

AICc(Phloro_mod1, Phloro_mod2, Phloro_mod3)
##             df     AICc
## Phloro_mod1  4 831.2629
## Phloro_mod2  5 832.2302
## Phloro_mod3  6 828.0473
summary(Phloro_mod3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: Phlorotannins ~ poly(Latitude, 3) + (1 | Site)
##    Data: Phlorotannins_metadata
## 
##      AIC      BIC   logLik deviance df.resid 
##    827.3    844.0   -407.7    815.3      114 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5999 -0.5788 -0.0462  0.2840  3.8222 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 0.06263  0.2503  
##  Residual             0.29291  0.5412  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value Pr(>|z|)    
## (Intercept)          2.6600     0.1165  22.838  < 2e-16 ***
## poly(Latitude, 3)1  -3.1963     1.2765  -2.504  0.01228 *  
## poly(Latitude, 3)2  -2.0276     1.2753  -1.590  0.11185    
## poly(Latitude, 3)3   3.8707     1.2767   3.032  0.00243 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1  0.000                
## ply(Ltt,3)2  0.001  0.000         
## ply(Ltt,3)3  0.000 -0.001  -0.013
r2(Phloro_mod3)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.546
##      Marginal R2: 0.435
predictions <- ggpredict(Phloro_mod3, terms ="Latitude [all]")
# Plot the GAM model predictions
Phloro_plot <- ggplot() +
  geom_jitter(data = Phlorotannins_metadata, aes(x= Latitude, y=Phlorotannins), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Phlorotannins (mg"~g^-1*" DW)", x="Latitude (°S)")  

Final plots

  1. Figure 2 - Abiotic predictors
Temperature_plot + Salinity_plot + Oxygen_plot + Phloro_plot+ plot_layout(ncol=2) 

Figure 2: Latitudinal variation of abiotic variables: temperature (°C), salinity (ppt) and dissolved oxygen (mgL) and visibility (m).

Save the plot in tiff format

ggsave("Latitudevsabiotic.jpeg", units="mm", width=125, height=85, dpi=300)
  1. Figure 3
Visibility_plot + Bare_plot+ Kelp_plot + UA_plot + Inv_plot + PredFish_plot + Size_plot + Richness_plot + PCOA_plot+ Herb_plot + plot_layout(ncol=2) 
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
## not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

Figure 3: Latitudinal variation of biotic and habitat variables

ggsave("FigurelatvsBio.jpeg", units="mm", width=110, height=140, dpi=300)

Part 4: Latitudinal patterns in consumption

Latitudinal patterns in consumption (predation and herbivory) were explored and analysed

4a): Predation vs latitude

First change latitude to positive numbers for plotting

Predation_env$Latitude <- Predation_env$Latitude*-1
C_24 <- filter(Predation_env, Bait2 == "Crabs" & Assay_time2 == 24)
C_2 <- filter(Predation_env, Bait2 == "Crabs" & Assay_time2 == 2)
S_24 <- filter(Predation_env, Bait2 == "Squidpops" & Assay_time2 == 24)
S_2 <- filter(Predation_env, Bait2 == "Squidpops" & Assay_time2 == 2)

1) Crabs 2h

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(C_2$Latitude, C_2$Consumed_individuals, xlab = "Latitude", ylab = "crabs consumed", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(C_2$Latitude, C_2$Consumed_individuals, pch = 21, col = mygrey, bg = mygrey)

descdist(C_2$Consumed_individuals, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  5 
## median:  1 
## mean:  1.958333 
## estimated sd:  2.043213 
## estimated skewness:  0.4960334 
## estimated kurtosis:  1.571122
normal_dist <- fitdist(C_2$Consumed_individuals, "exp")
plot(normal_dist)

C2_mod1 <- glmer(Consumed_individuals ~ Latitude+ (1|Site),  data= C_2, family = poisson())
C2_mod2 <- glmer(Consumed_individuals ~poly(Latitude,2)+ (1|Site),  data= C_2, family = poisson())
C2_mod3 <- glmer(Consumed_individuals ~poly(Latitude,3)+ (1|Site),  data= C_2, family = poisson())
check_model(C2_mod1)

AICc(C2_mod2, C2_mod1, C2_mod3)
##         df     AICc
## C2_mod2  4 424.9806
## C2_mod1  3 423.1983
## C2_mod3  5 425.7216
BIC(C2_mod2, C2_mod1, C2_mod3)
##         df      BIC
## C2_mod2  4 435.7828
## C2_mod1  3 431.3539
## C2_mod3  5 439.1328
summary(C2_mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Latitude + (1 | Site)
##    Data: C_2
## 
##      AIC      BIC   logLik deviance df.resid 
##    423.0    431.4   -208.5    417.0      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92403 -0.94180  0.03414  0.62864  3.02889 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.414    0.6434  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  2.33581    0.91197   2.561   0.0104 *
## Latitude    -0.06553    0.03032  -2.161   0.0307 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## Latitude -0.975
r2(C2_mod1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.601
##      Marginal R2: 0.198

2) Crabs 24h

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(C_24$Latitude, C_24$Consumed_individuals, xlab = "Latitude", ylab = "crabs consumed 24 h", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(C_24$Latitude, C_24$Consumed_individuals, pch = 21, col = mygrey, bg = mygrey)

descdist(C_24$Consumed_individuals, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  5 
## median:  5 
## mean:  3.791667 
## estimated sd:  1.511279 
## estimated skewness:  -0.8418895 
## estimated kurtosis:  2.211133
normal_dist <- fitdist(C_24$Consumed_individuals, "exp")
plot(normal_dist)

C24_mod1 <- glm(Consumed_individuals ~ Latitude,  data= C_24, family = poisson())
C24_mod2 <- glm(Consumed_individuals ~poly(Latitude,2),  data= C_24, family = poisson())
C24_mod3 <- glm(Consumed_individuals ~poly(Latitude,3),  data= C_24, family = poisson())
check_overdispersion(C24_mod1)
## # Overdispersion test
## 
##        dispersion ratio =  0.420
##   Pearson's Chi-Squared = 49.526
##                 p-value =      1
## No overdispersion detected.
check_model(C24_mod1)

AICc(C24_mod1, C24_mod2, C24_mod3)
##          df     AICc
## C24_mod1  2 433.7401
## C24_mod2  3 434.2198
## C24_mod3  4 436.0433
BIC(C24_mod1, C24_mod2, C24_mod3)
##          df      BIC
## C24_mod1  2 439.2125
## C24_mod2  3 442.3754
## C24_mod3  4 446.8454
summary(C24_mod1)
## 
## Call:
## glm(formula = Consumed_individuals ~ Latitude, family = poisson(), 
##     data = C_24)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.452269   0.208192  11.779  < 2e-16 ***
## Latitude    -0.038689   0.007222  -5.357 8.46e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 58.265  on 118  degrees of freedom
## AIC: 433.64
## 
## Number of Fisher Scoring iterations: 4
r2(C24_mod1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.424

Calculate predictions for Crab consumption after 2 and 24 h

predictions <- ggpredict(C2_mod1, terms ="Latitude")
predictions1 <- ggpredict(C24_mod1, terms ="Latitude")

Add new column to mix datasets

predictions$Assay_time <- "2 h"
predictions1$Assay_time <- "24 h"
C_2$Assay_time <- "2 h"
C_24$Assay_time <- "24 h"

#Bind datasets

total <- rbind(predictions, predictions1)

total2 <- rbind (C_2, C_24)

3) Crabs 2 and 24 h prediction plot

ggplot()+
  geom_ribbon(data = total, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, group = Assay_time, fill = Assay_time ), alpha = .25)  +
  geom_line(data = total, aes(x = x, y = predicted, color = Assay_time, group = Assay_time), size = 1) +
  geom_jitter(data = total2, aes(x= Latitude, y=Consumed_individuals, color =Assay_time), width=0.15, height = 0.2, size = 1, alpha = 0.5)+
  labs(x = "Latitude", y = "Crabs consumed")  + Theme1 +
  scale_color_manual("Assay time",values=c("darkslateblue", "orange1"))+
   scale_fill_manual("Assay time", values=c("darkslateblue", "orange1")) + 
  theme(legend.position = c(0.86, 0.81),legend.title = element_text(size = 8), 
               legend.text = element_text(size = 6.3), legend.background = element_rect(fill="NA",
                                         size=1, 
                                         color ="NA")) +
  scale_x_continuous(breaks = seq(15, 40, 5)) +
  scale_y_continuous(limits = c(0,6.5), breaks = seq(0, 5, 1))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 27 rows containing missing values or values outside the scale range
## (`geom_point()`).

Save plot in tiff format

ggsave("Crabs.jpeg", units="mm", width=100, height=65, dpi=300)
## Warning: Removed 19 rows containing missing values or values outside the scale range
## (`geom_point()`).

4) Squidpops 2h

a. Analyses

descdist(S_2$Consumed_individuals, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  5 
## median:  2 
## mean:  2.375 
## estimated sd:  1.974895 
## estimated skewness:  0.1159755 
## estimated kurtosis:  1.430103
normal_dist <- fitdist(S_2$Consumed_individuals, "nbinom")
plot(normal_dist)

S2_mod1 <- glmer(Consumed_individuals ~ Latitude + (1|Site),  data= S_2, family = poisson())
S2_mod2 <- glmer(Consumed_individuals ~poly(Latitude,2)+ (1|Site),  data= S_2, family = poisson())
S2_mod3 <- glmer(Consumed_individuals ~poly(Latitude,3)+ (1|Site),  data= S_2, family = poisson())
check_overdispersion(S2_mod1)
## # Overdispersion test
## 
##        dispersion ratio =   1.050
##   Pearson's Chi-Squared = 122.805
##                 p-value =   0.338
## No overdispersion detected.
check_model(S2_mod1)

AICc(S2_mod2, S2_mod1, S2_mod3)
##         df     AICc
## S2_mod2  4 462.0740
## S2_mod1  3 460.3905
## S2_mod3  5 463.5723
BIC(S2_mod2, S2_mod1, S2_mod3)
##         df      BIC
## S2_mod2  4 472.8761
## S2_mod1  3 468.5461
## S2_mod3  5 476.9835
summary(S2_mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Latitude + (1 | Site)
##    Data: S_2
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.2    468.5   -227.1    454.2      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7204 -0.8067 -0.0802  0.7098  2.3659 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5241   0.7239  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.14773    0.99008   0.149    0.881
## Latitude     0.01733    0.03209   0.540    0.589
## 
## Correlation of Fixed Effects:
##          (Intr)
## Latitude -0.974
r2(S2_mod1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.618
##      Marginal R2: 0.016

5) Squidpops 24h

a. Analyses

descdist(S_24$Consumed_individuals, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  5 
## median:  5 
## mean:  4.25 
## estimated sd:  1.258584 
## estimated skewness:  -1.822846 
## estimated kurtosis:  5.495426
normal_dist <- fitdist(S_24$Consumed_individuals, "pois")
plot(normal_dist)

S24_mod1 <- glm(Consumed_individuals ~ Latitude,  data= S_24, family = poisson())
S24_mod2 <- glm(Consumed_individuals ~poly(Latitude,2),  data= S_24, family = poisson())
S24_mod3 <- glm(Consumed_individuals ~poly(Latitude,3),  data= S_24, family = poisson())
check_overdispersion(S24_mod1)
## # Overdispersion test
## 
##        dispersion ratio =  0.348
##   Pearson's Chi-Squared = 41.044
##                 p-value =      1
## No overdispersion detected.
check_model(S24_mod1)

AICc(S24_mod2, S24_mod1, S24_mod3)
##          df     AICc
## S24_mod2  3 447.7243
## S24_mod1  2 446.1367
## S24_mod3  4 449.8652
BIC(S24_mod2, S24_mod1, S24_mod3)
##          df      BIC
## S24_mod2  3 455.8799
## S24_mod1  2 451.6091
## S24_mod3  4 460.6673
summary(S24_mod1)
## 
## Call:
## glm(formula = Consumed_individuals ~ Latitude, family = poisson(), 
##     data = S_24)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.972899   0.196634  10.033  < 2e-16 ***
## Latitude    -0.017890   0.006608  -2.707  0.00678 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 54.280  on 118  degrees of freedom
## AIC: 446.03
## 
## Number of Fisher Scoring iterations: 4
r2(S24_mod1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.149
#For squidpops 24 h 
predictions3 <- ggpredict(S24_mod1, terms ="Latitude")

Add new column to mix datasets

S_2$Assay_time <- "2 h"
S_24$Assay_time <- "24 h"
predictions3$Assay_time <- "24 h"
##Keep only colums for plotting
total2 <- rbind (S_2, S_24)

3) Squidpops 2 and 24 h prediction plot

ggplot() +
  geom_ribbon(data = predictions3, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, group = Assay_time, fill = Assay_time ), alpha = .25) +
  geom_line(data = predictions3, aes(x = x, y = predicted, color = Assay_time, group = Assay_time ), size = 1) +
  geom_jitter(data = total2, aes(x= Latitude, y=Consumed_individuals, color =Assay_time), width=0.15, height = 0.2, size = 1, alpha = .5)+
  labs(x = "Latitude", y = "Squidpops consumed")  + Theme1 +
  scale_color_manual("Assay time",values=c("darkslateblue", "orange1"))+
   scale_fill_manual("Assay time", values=c("orange1")) + 
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(15, 40, 5)) +
  scale_y_continuous(limits = c(0,6.5), breaks = seq(0, 5, 1))
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

Save plot in tiff format

ggsave("Squidpops.jpeg", units="mm", width=100, height=65, dpi=300)
## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).

4b): Herbivory vs latitude

Herbivory_env$Latitude <- Herbivory_env$Latitude*-1

a. Analyses

descdist(Herbivory_env$Consumed_proportion, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  1 
## median:  0 
## mean:  0.2260169 
## estimated sd:  0.3698909 
## estimated skewness:  1.239303 
## estimated kurtosis:  2.780775
normal_dist <- fitdist(Herbivory_env$Consumed_proportion, "exp")
plot(normal_dist)

Herb_mod1 <- glmmTMB(Consumed_proportion ~ Latitude+ (1|Site),  data= Herbivory_env, family = tweedie())
Herb_mod2 <- glmmTMB(Consumed_proportion ~poly(Latitude,2)+ (1|Site),  data= Herbivory_env, family = tweedie())
Herb_mod3 <- glmmTMB(Consumed_proportion ~poly(Latitude,3)+ (1|Site),  data= Herbivory_env, family = tweedie())
check_model(Herb_mod1)
## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(Herb_mod2, Herb_mod1, Herb_mod3)
##           df     AICc
## Herb_mod2  6 732.3190
## Herb_mod1  5 730.7644
## Herb_mod3  7 732.0831
summary(Herb_mod1)
##  Family: tweedie  ( log )
## Formula:          Consumed_proportion ~ Latitude + (1 | Site)
## Data: Herbivory_env
## 
##      AIC      BIC   logLik deviance df.resid 
##    730.7    752.6   -360.3    720.7      585 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.69     2.586   
## Number of obs: 590, groups:  Site, 12
## 
## Dispersion parameter for tweedie family (): 0.875 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -8.9565     3.5611  -2.515   0.0119 *
## Latitude      0.1899     0.1138   1.668   0.0953 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(Herb_mod1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.956
##      Marginal R2: 0.195

b. Plot predicted relationship

ggplot() +
  geom_jitter(data = Herbivory_env, aes(x= Latitude, y=Consumed_proportion), color = "orange1", width=0.2, height = 0.02, size = 1, alpha = .5)+
  labs(x = "Latitude (°S)", y = "Blade consumed proportion")  + Theme1  + 
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(15, 40, 5)) +
  scale_y_continuous(limits = c(-0.1,1.2), breaks = seq(0, 1, 0.2))

Save plot to tiff format

ggsave("Seaweed24.jpeg", units="mm", width=105, height=65, dpi=300)

Part 5: Predation intensity: Effect of site-specific variables

Here, the individual effects of each predictor will be tested.

  1. First, create datasets to be used in the analysis
Predation_analysis <- Predation_6 %>%
  dplyr::select(-Rope, -Bait, -Camera, -Assay_time, -Initial_Individuals, -Remaining_individuals, -Survival_proportion,-Consumed_proportion)

#Separate by Bait and Assay time 

Crabs_24 <- filter(Predation_analysis, Bait2 == "Crabs" & Assay_time2 == 24)
Crabs_2 <- filter(Predation_analysis, Bait2 == "Crabs" & Assay_time2 == 2)
Squidpops_24 <- filter(Predation_analysis, Bait2 == "Squidpops" & Assay_time2 == 24)
Squidpops_2 <- filter(Predation_analysis, Bait2 == "Squidpops" & Assay_time2 == 2)

5a): Standardisation and test autocorrelation

1) Standardise

#Crabs 24h
Crabs24_num <- Crabs_24 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude, -Day_length) 
Crabs24_stand <- decostand(Crabs24_num, method="range", na.rm = TRUE)
#Test autocorrelation 
corrplot(cor(Crabs24_stand), method = "number")

No multicolinearity found

Crabs24_stand$Site <- Crabs_24$Site
Crabs24_stand$Consumed_individuals<- Crabs_24$Consumed_individuals
Crabs24_stand$Latitude<- Crabs_24$Latitude
#Crabs 2h 
Crabs2_num <- Crabs_2 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude) 
Crabs2_stand <- decostand(Crabs2_num, method="range", na.rm = TRUE)
Crabs2_stand$Site <- Crabs_2$Site
Crabs2_stand$Consumed_individuals<- Crabs_2$Consumed_individuals
Crabs2_stand$Latitude<- Crabs_2$Latitude
#Squidpops 24h
Squidpops24_num <- Squidpops_24 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude) 
Squidpops24_stand <- decostand(Squidpops24_num, method="range", na.rm = TRUE)
Squidpops24_stand$Site <- Squidpops_24$Site
Squidpops24_stand$Consumed_individuals<- Squidpops_24$Consumed_individuals
Squidpops24_stand$Latitude<- Squidpops_24$Latitude
#Squidpops 2h
Squidpops2_num <- Squidpops_2 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude) 
Squidpops2_stand <- decostand(Squidpops2_num, method="range", na.rm = TRUE)
Squidpops2_stand$Site <- Squidpops_2$Site
Squidpops2_stand$Consumed_individuals<- Squidpops_2$Consumed_individuals
Squidpops2_stand$Latitude<- Squidpops_2$Latitude

5b) Model selection: Crabs 2h

  1. Temperature
my_glmT1<- glmer(Consumed_individuals~  Temperature + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmT1)

AICc(my_glmT1)
## [1] 420.2918
summary(my_glmT1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    420.1    428.4   -207.0    414.1      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93381 -0.92312  0.06015  0.65176  2.79744 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.2979   0.5458  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -0.8086     0.4487  -1.802  0.07151 . 
## Temperature   1.8546     0.6104   3.038  0.00238 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temperature -0.918
r2(my_glmT1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.599
##      Marginal R2: 0.307
  1. Salinity
my_glmSC2_1<- glmer(Consumed_individuals~  Salinity_ppt + (1|Site), data=Crabs2_stand, family = poisson)


check_overdispersion(my_glmSC2_1)
## # Overdispersion test
## 
##        dispersion ratio =   1.189
##   Pearson's Chi-Squared = 139.141
##                 p-value =    0.08
## No overdispersion detected.
check_model(my_glmSC2_1)

AICc(my_glmSC2_1)
## [1] 424.5386
summary(my_glmSC2_1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Salinity_ppt + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    424.3    432.7   -209.2    418.3      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92849 -0.91521  0.04145  0.66421  2.99635 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.468    0.6841  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -0.6810     0.6667  -1.022   0.3070  
## Salinity_ppt   1.5821     0.9192   1.721   0.0852 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Salinty_ppt -0.947
r2(my_glmSC2_1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.603
##      Marginal R2: 0.148
  1. Dissolved oxygen
my_glmOxC2_1<- glmer(Consumed_individuals~  D_oxygen_mgL + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmOxC2_1)

AICc(my_glmOxC2_1)
## [1] 422.5227
summary(my_glmOxC2_1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ D_oxygen_mgL + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    422.3    430.7   -208.2    416.3      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9665 -0.9020  0.0611  0.5761  3.2095 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.365    0.6042  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -0.2934     0.3526  -0.832   0.4054  
## D_oxygen_mgL   1.5663     0.6451   2.428   0.0152 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## D_oxygn_mgL -0.836
r2(my_glmOxC2_1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.585
##      Marginal R2: 0.215
  1. Visibility
my_glmVC2_1<- glmer(Consumed_individuals~V_Visibility + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmVC2_1)

AICc(my_glmVC2_1)
## [1] 427.2626
summary(my_glmVC2_1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ V_Visibility + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    427.1    435.4   -210.5    421.1      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93460 -0.89601  0.07717  0.53694  3.02786 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5924   0.7697  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.39200    0.44714   0.877    0.381
## V_Visibility -0.01098    0.89880  -0.012    0.990
## 
## Correlation of Fixed Effects:
##             (Intr)
## V_Visibilty -0.847
r2(my_glmVC2_1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.592
##      Marginal R2: 0.000
  1. Community composition
my_glmC2_PCO1<- glmer(Consumed_individuals~PCo1 + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_PCO1)

AICc(my_glmC2_PCO1)
## [1] 426.871
summary(my_glmC2_PCO1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ PCo1 + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    426.7    435.0   -210.3    420.7      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92572 -0.90531  0.06451  0.54796  2.99751 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5818   0.7628  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.5359     0.3320   1.614    0.107
## PCo1         -0.4240     0.6764  -0.627    0.531
## 
## Correlation of Fixed Effects:
##      (Intr)
## PCo1 -0.703
r2(my_glmC2_PCO1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.597
##      Marginal R2: 0.023
  1. Kelp Cover
my_glmC2_K1<- glmer(Consumed_individuals~Kelp + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_K1)

AICc(my_glmC2_K1)
## [1] 424.3606
summary(my_glmC2_K1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Kelp + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    424.2    432.5   -209.1    418.2      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92113 -0.89423  0.02569  0.55099  3.16011 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4494   0.6704  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.02724    0.29757   0.092   0.9271  
## Kelp         1.32373    0.73207   1.808   0.0706 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## Kelp -0.705
r2(my_glmC2_K1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.589
##      Marginal R2: 0.138
  1. Bare rock cover
my_glmC2_BR1<- glmer(Consumed_individuals~Bare_rock + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_BR1)

AICc(my_glmC2_BR1)
## [1] 424.1326
summary(my_glmC2_BR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Bare_rock + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    423.9    432.3   -209.0    417.9      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9360 -0.9040  0.0409  0.5802  3.2149 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4321   0.6573  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.8447     0.3094   2.730  0.00633 **
## Bare_rock    -1.1994     0.6315  -1.899  0.05751 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Bare_rock -0.741
r2(my_glmC2_BR1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.593
##      Marginal R2: 0.163
  1. Understorey algae cover
my_glmC2_UA1<- glmer(Consumed_individuals~Understorey_algae + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_UA1)

AICc(my_glmC2_UA1)
## [1] 427.0992
summary(my_glmC2_UA1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Understorey_algae + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    426.9    435.3   -210.4    420.9      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92743 -0.90039  0.08966  0.54638  3.04101 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5829   0.7635  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)         0.4574     0.2908   1.573    0.116
## Understorey_algae  -0.2698     0.6637  -0.406    0.684
## 
## Correlation of Fixed Effects:
##             (Intr)
## Undrstry_lg -0.584
r2(my_glmC2_UA1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.592
##      Marginal R2: 0.009
  1. Species richness
my_glmC2_SR1<- glmer(Consumed_individuals~Mean_SR + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_SR1)

AICc(my_glmC2_SR1)
## [1] 425.1803
summary(my_glmC2_SR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Mean_SR + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    425.0    433.3   -209.5    419.0      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.90607 -0.90370  0.06207  0.59243  3.08960 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4998   0.707   
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.1651     0.4393  -0.376    0.707
## Mean_SR       1.1522     0.7741   1.488    0.137
## 
## Correlation of Fixed Effects:
##         (Intr)
## Mean_SR -0.864
r2(my_glmC2_SR1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.598
##      Marginal R2: 0.105
  1. Invertebrate predators
my_glmC2_Inv1<- glmer(Consumed_individuals~Invertebrate_predators + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_Inv1)

AICc(my_glmC2_Inv1)
## [1] 427.1083
summary(my_glmC2_Inv1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Invertebrate_predators + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    426.9    435.3   -210.5    420.9      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93478 -0.89816  0.08557  0.54475  3.04770 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5865   0.7658  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)
## (Intercept)              0.2603     0.4017   0.648    0.517
## Invertebrate_predators   0.2879     0.7302   0.394    0.693
## 
## Correlation of Fixed Effects:
##             (Intr)
## Invrtbrt_pr -0.808
r2(my_glmC2_Inv1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.593
##      Marginal R2: 0.009
  1. Predatory fish abundance
my_glmC2_PF1<- glmer(Consumed_individuals~Predatory_fish + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_PF1)

AICc(my_glmC2_PF1)
## [1] 426.858
summary(my_glmC2_PF1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Predatory_fish + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    426.7    435.0   -210.3    420.7      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92759 -0.89488  0.07857  0.55013  3.05066 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5723   0.7565  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)      0.2205     0.3528   0.625    0.532
## Predatory_fish   0.4397     0.6854   0.642    0.521
## 
## Correlation of Fixed Effects:
##             (Intr)
## Predtry_fsh -0.748
r2(my_glmC2_PF1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.593
##      Marginal R2: 0.022
  1. Fish size
my_glmC2_FL1<- glmer(Consumed_individuals~Fish_size + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_FL1)

AICc(my_glmC2_FL1)
## [1] 426.9261
summary(my_glmC2_FL1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Fish_size + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    426.7    435.1   -210.4    420.7      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93714 -0.91055  0.06116  0.54889  3.03983 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5781   0.7603  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.06379    0.60517   0.105    0.916
## Fish_size    0.52233    0.89451   0.584    0.559
## 
## Correlation of Fixed Effects:
##           (Intr)
## Fish_size -0.921
r2(my_glmC2_FL1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.594
##      Marginal R2: 0.020

5c) Model selection: Crabs 24h

  1. Temperature
my_glmC24_T1<- glmer(Consumed_individuals~Temperature + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_T1)

AICc(my_glmC24_T1)
## [1] 438.9369
summary(my_glmC24_T1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    438.7    447.1   -216.4    432.7      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1336 -0.2824  0.1619  0.4568  1.3189 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev.
##  Site   (Intercept) 0.0003175 0.01782 
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.7379     0.1353   5.456 4.88e-08 ***
## Temperature   0.8747     0.1776   4.925 8.42e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temperature -0.936
r2(my_glmC24_T1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.230
##      Marginal R2: 0.229
  1. Salinity
my_glmC24_S1<- glmer(Consumed_individuals~Salinity_ppt + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_S1)

AICc(my_glmC24_S1)
## [1] 443.1697
summary(my_glmC24_S1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Salinity_ppt + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    443.0    451.3   -218.5    437.0      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.01001 -0.46529  0.09323  0.42898  1.47247 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.01181  0.1087  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.6358     0.1945   3.269 0.001078 ** 
## Salinity_ppt   0.9897     0.2621   3.777 0.000159 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Salinty_ppt -0.956
r2(my_glmC24_S1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.233
##      Marginal R2: 0.195
  1. Dissolved oxygen
my_glmC24_Ox1<- glmer(Consumed_individuals~D_oxygen_mgL + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_Ox1)

AICc(my_glmC24_Ox1)
## [1] 452.5543
summary(my_glmC24_Ox1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ D_oxygen_mgL + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    452.3    460.7   -223.2    446.3      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0829 -0.3910  0.1934  0.3175  1.2599 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.05924  0.2434  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.1962     0.1534   7.796 6.38e-15 ***
## D_oxygen_mgL   0.2413     0.2866   0.842      0.4    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## D_oxygn_mgL -0.831
r2(my_glmC24_Ox1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.216
##      Marginal R2: 0.017
  1. Visibility
my_glmC24_V1<- glmer(Consumed_individuals~V_Visibility + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_V1)

AICc(my_glmC24_V1)
## [1] 452.5114
summary(my_glmC24_V1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ V_Visibility + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    452.3    460.7   -223.2    446.3      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0495 -0.3868  0.2021  0.3902  1.2724 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.05948  0.2439  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.4210     0.1604   8.858   <2e-16 ***
## V_Visibility  -0.2825     0.3267  -0.865    0.387    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## V_Visibilty -0.846
r2(my_glmC24_V1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.218
##      Marginal R2: 0.019
  1. Consumer community
my_glmC24_PCO1<- glmer(Consumed_individuals~PCo1 + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_PCO1)

AICc(my_glmC24_PCO1)
## [1] 450.5558
summary(my_glmC24_PCO1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ PCo1 + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    450.3    458.7   -222.2    444.3      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9758 -0.3167  0.2087  0.3774  1.2274 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.04741  0.2177  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.4411     0.1109  12.993   <2e-16 ***
## PCo1         -0.3920     0.2300  -1.705   0.0883 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## PCo1 -0.698
r2(my_glmC24_PCO1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.223
##      Marginal R2: 0.065
  1. Kelp cover
my_glmC24_K1<- glmer(Consumed_individuals~Kelp + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_K1)

AICc(my_glmC24_K1)
## [1] 448.7633
summary(my_glmC24_K1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Kelp + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    448.6    456.9   -221.3    442.6      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0196 -0.4200  0.1414  0.4214  1.2459 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.03546  0.1883  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1427     0.1040  10.987   <2e-16 ***
## Kelp          0.5844     0.2518   2.321   0.0203 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## Kelp -0.715
r2(my_glmC24_K1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.210
##      Marginal R2: 0.090
  1. Bare rock cover
my_glmC24_BR1<- glmer(Consumed_individuals~Bare_rock + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_BR1)

AICc(my_glmC24_BR1)
## [1] 447.5098
summary(my_glmC24_BR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Bare_rock + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    447.3    455.7   -220.7    441.3      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0745 -0.3637  0.1577  0.3357  1.5700 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.02761  0.1662  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.51816    0.09993  15.192  < 2e-16 ***
## Bare_rock   -0.56632    0.20938  -2.705  0.00683 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Bare_rock -0.733
r2(my_glmC24_BR1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.216
##      Marginal R2: 0.123
  1. Understorey algae cover
my_glmC24_UA1<- glmer(Consumed_individuals~Understorey_algae + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_UA1)

AICc(my_glmC24_UA1)
## [1] 453.1994
summary(my_glmC24_UA1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Understorey_algae + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    453.0    461.4   -223.5    447.0      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0561 -0.4182  0.1760  0.3756  1.3128 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.06408  0.2531  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        1.28831    0.10876  11.846   <2e-16 ***
## Understorey_algae  0.05185    0.24524   0.211    0.833    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Undrstry_lg -0.590
r2(my_glmC24_UA1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.216
##      Marginal R2: 0.001
  1. Species richness
my_glmC24_SR1<- glmer(Consumed_individuals~Mean_SR + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_SR1)

AICc(my_glmC24_SR1)
## [1] 451.6355
summary(my_glmC24_SR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Mean_SR + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    451.4    459.8   -222.7    445.4      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9950 -0.3586  0.1690  0.4479  1.2775 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.05277  0.2297  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1192     0.1647   6.797 1.07e-11 ***
## Mean_SR       0.3830     0.2925   1.309     0.19    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## Mean_SR -0.867
r2(my_glmC24_SR1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.216
##      Marginal R2: 0.040
  1. Abundance of invertebrate predators
my_glmC24_Inv1<- glmer(Consumed_individuals~Invertebrate_predators + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_Inv1)

AICc(my_glmC24_Inv1)
## [1] 453.0005
summary(my_glmC24_Inv1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Invertebrate_predators + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    452.8    461.2   -223.4    446.8      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0448 -0.4524  0.1776  0.4005  1.3483 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.06186  0.2487  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              1.3610     0.1462   9.308   <2e-16 ***
## Invertebrate_predators  -0.1345     0.2703  -0.498    0.619    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Invrtbrt_pr -0.805
r2(my_glmC24_Inv1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.214
##      Marginal R2: 0.007
  1. Abundance of fish predators
my_glmC24_PF1<- glmer(Consumed_individuals~Predatory_fish + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_PF1)

AICc(my_glmC24_PF1)
## [1] 452.9183
summary(my_glmC24_PF1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Predatory_fish + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    452.7    461.1   -223.4    446.7      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0344 -0.3964  0.1898  0.4233  1.2267 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.06197  0.2489  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      1.2463     0.1308   9.528   <2e-16 ***
## Predatory_fish   0.1465     0.2549   0.575    0.565    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Predtry_fsh -0.748
r2(my_glmC24_PF1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.216
##      Marginal R2: 0.008
  1. Fish legth (cm)
my_glmC24_FL1<- glmer(Consumed_individuals~Fish_size + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_FL1)

AICc(my_glmC24_FL1)
## [1] 448.4084
summary(my_glmC24_FL1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Fish_size + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    448.2    456.6   -221.1    442.2      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0671 -0.2882  0.2492  0.3076  1.0942 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.03566  0.1888  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.8921     0.1932   4.618 3.88e-06 ***
## Fish_size     0.6622     0.2796   2.369   0.0179 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Fish_size -0.926
r2(my_glmC24_FL1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.226
##      Marginal R2: 0.108

5e): Model selection: Squidpops 2h

  1. Temperature
my_glmS2_T1<- glmer(Consumed_individuals~Temperature + (1|Site), data=Squidpops2_stand, family = poisson)

check_overdispersion(my_glmS2_T1)
## # Overdispersion test
## 
##        dispersion ratio =   1.054
##   Pearson's Chi-Squared = 123.348
##                 p-value =   0.326
## No overdispersion detected.
check_model(my_glmS2_T1)

AICc(my_glmS2_T1)
## [1] 460.4034
summary(my_glmS2_T1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.2    468.6   -227.1    454.2      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7049 -0.8189 -0.1183  0.6789  2.3417 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5106   0.7145  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.4209     0.5139   0.819    0.413
## Temperature   0.3832     0.7173   0.534    0.593
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temperature -0.903
r2(my_glmS2_T1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.612
##      Marginal R2: 0.016
  1. Salinity
my_glmS2_S1<- glmer(Consumed_individuals~Salinity_ppt + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_S1)

AICc(my_glmS2_S1)
## [1] 460.4596
summary(my_glmS2_S1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Salinity_ppt + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.3    468.6   -227.1    454.3      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.72415 -0.80815 -0.08599  0.69023  2.37161 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5202   0.7212  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)    0.9474     0.6296   1.505    0.132
## Salinity_ppt  -0.4192     0.8834  -0.475    0.635
## 
## Correlation of Fixed Effects:
##             (Intr)
## Salinty_ppt -0.935
r2(my_glmS2_S1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.615
##      Marginal R2: 0.012
  1. Dissolved oxygen
my_glmS2_Ox1<- glmer(Consumed_individuals~D_oxygen_mgL + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_Ox1)

AICc(my_glmS2_Ox1)
## [1] 457.1796
summary(my_glmS2_Ox1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ D_oxygen_mgL + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    457.0    465.3   -225.5    451.0      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.72360 -0.79117 -0.07451  0.73700  2.40980 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.3718   0.6098  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    0.1137     0.3581   0.318   0.7508  
## D_oxygen_mgL   1.2713     0.6529   1.947   0.0515 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## D_oxygn_mgL -0.844
r2(my_glmS2_Ox1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.606
##      Marginal R2: 0.166
  1. Visibility
my_glmS2_V1<- glmer(Consumed_individuals~V_Visibility + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_V1)

AICc(my_glmS2_V1)
## [1] 460.1753
summary(my_glmS2_V1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ V_Visibility + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.0    468.3   -227.0    454.0      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7021 -0.8107 -0.1061  0.7043  2.4274 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5109   0.7148  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)    0.4194     0.4157   1.009    0.313
## V_Visibility   0.5814     0.8141   0.714    0.475
## 
## Correlation of Fixed Effects:
##             (Intr)
## V_Visibilty -0.847
r2(my_glmS2_V1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.617
##      Marginal R2: 0.028
  1. Community composition
my_glmS2_PCO1<- glmer(Consumed_individuals~PCo1 + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_PCO1)

AICc(my_glmS2_PCO1)
## [1] 457.9922
summary(my_glmS2_PCO1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ PCo1 + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    457.8    466.1   -225.9    451.8      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.74367 -0.78077 -0.08552  0.71274  2.46497 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.3944   0.628   
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.3370     0.2866   1.176   0.2396  
## PCo1          0.9343     0.5457   1.712   0.0869 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## PCo1 -0.726
r2(my_glmS2_PCO1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.604
##      Marginal R2: 0.133
  1. Kelp cover
my_glmS2_K1<- glmer(Consumed_individuals~Kelp + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_K1)

AICc(my_glmS2_K1)
## [1] 459.4048
summary(my_glmS2_K1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Kelp + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    459.2    467.6   -226.6    453.2      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6926 -0.8389 -0.1096  0.6556  2.4268 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4541   0.6739  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.4384     0.2948   1.487    0.137
## Kelp          0.8480     0.7294   1.163    0.245
## 
## Correlation of Fixed Effects:
##      (Intr)
## Kelp -0.702
r2(my_glmS2_K1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.606
##      Marginal R2: 0.067
  1. Bare rock cover
my_glmS2_BR1<- glmer(Consumed_individuals~Bare_rock + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_BR1)

AICc(my_glmS2_BR1)
## [1] 460.4304
summary(my_glmS2_BR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Bare_rock + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.2    468.6   -227.1    454.2      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7003 -0.8230 -0.1211  0.6796  2.3980 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5128   0.7161  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.7914     0.3284   2.410    0.016 *
## Bare_rock    -0.3255     0.6415  -0.507    0.612  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Bare_rock -0.738
r2(my_glmS2_BR1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.612
##      Marginal R2: 0.014
  1. Understorey algae cover
my_glmS2_UA1<- glmer(Consumed_individuals~Understorey_algae + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_UA1)

AICc(my_glmS2_UA1)
## [1] 460.6711
summary(my_glmS2_UA1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Understorey_algae + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.5    468.8   -227.2    454.5      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.71045 -0.81805 -0.09991  0.68648  2.36196 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5282   0.7267  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        0.68441    0.27845   2.458    0.014 *
## Understorey_algae -0.06958    0.62237  -0.112    0.911  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Undrstry_lg -0.592
r2(my_glmS2_UA1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.614
##      Marginal R2: 0.001
  1. Species richness
my_glmS2_SR1<- glmer(Consumed_individuals~Mean_SR + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_SR1)

AICc(my_glmS2_SR1)
## [1] 460.6386
summary(my_glmS2_SR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Mean_SR + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.4    468.8   -227.2    454.4      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.71377 -0.81473 -0.09892  0.68827  2.36015 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5293   0.7275  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.7440     0.4298   1.731   0.0835 .
## Mean_SR      -0.1641     0.7736  -0.212   0.8320  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## Mean_SR -0.853
r2(my_glmS2_SR1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.615
##      Marginal R2: 0.003
  1. Abundance of invertebrate predators
my_glmS2_I1<- glmer(Consumed_individuals~Invertebrate_predators + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_I1)

AICc(my_glmS2_I1)
## [1] 460.0393
summary(my_glmS2_I1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Invertebrate_predators + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    459.8    468.2   -226.9    453.8      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7185 -0.8250 -0.1281  0.7024  2.4144 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4876   0.6983  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)  
## (Intercept)              0.9020     0.3548   2.542    0.011 *
## Invertebrate_predators  -0.5317     0.6507  -0.817    0.414  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Invrtbrt_pr -0.792
r2(my_glmS2_I1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.609
##      Marginal R2: 0.036
  1. Abundance of fish predators
my_glmS2_PF1<- glmer(Consumed_individuals~Predatory_fish + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_PF1)

AICc(my_glmS2_PF1)
## [1] 460.5443
summary(my_glmS2_PF1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Predatory_fish + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.3    468.7   -227.2    454.3      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.70479 -0.82195 -0.09594  0.67244  2.34183 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5245   0.7242  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)  
## (Intercept)      0.5755     0.3320   1.733    0.083 .
## Predatory_fish   0.2390     0.6392   0.374    0.709  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Predtry_fsh -0.739
r2(my_glmS2_PF1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.615
##      Marginal R2: 0.008
  1. Fish length
my_glmS2_FL1<- glmer(Consumed_individuals~Fish_size + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_FL1)

AICc(my_glmS2_FL1)
## [1] 458.2359
summary(my_glmS2_FL1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Fish_size + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    458.0    466.4   -226.0    452.0      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.72996 -0.80127 -0.09444  0.71113  2.33612 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4281   0.6543  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   1.4069     0.4953   2.841   0.0045 **
## Fish_size    -1.2026     0.7565  -1.590   0.1119   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Fish_size -0.911
r2(my_glmS2_FL1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.617
##      Marginal R2: 0.124

5f) Model selection: Squidpops 24h

  1. Temperature
my_glmS24_T1<- glmer(Consumed_individuals~Temperature+ (1|Site), data=Squidpops24_stand, family = poisson)
## boundary (singular) fit: see help('isSingular')
my_glmS24_T1_v2<- glm(Consumed_individuals~Temperature, data=Squidpops24_stand, family = poisson)
check_overdispersion(my_glmS24_T1_v2)
## # Overdispersion test
## 
##        dispersion ratio =  0.331
##   Pearson's Chi-Squared = 39.057
##                 p-value =      1
## No overdispersion detected.
check_model(my_glmS24_T1_v2)

AICc(my_glmS24_T1_v2)
## [1] 443.5292
summary(my_glmS24_T1_v2)
## 
## Call:
## glm(formula = Consumed_individuals ~ Temperature, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1227     0.1166   9.632  < 2e-16 ***
## Temperature   0.4879     0.1575   3.098  0.00195 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 51.673  on 118  degrees of freedom
## AIC: 443.43
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_T1_v2)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.200
  1. Salinity
my_glmS24_S1<- glm(Consumed_individuals~Salinity_ppt, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_S1)

AICc(my_glmS24_S1)
## [1] 447.5637
summary(my_glmS24_S1)
## 
## Call:
## glm(formula = Consumed_individuals ~ Salinity_ppt, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.1271     0.1434   7.857 3.93e-15 ***
## Salinity_ppt   0.4670     0.1956   2.387    0.017 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 55.707  on 118  degrees of freedom
## AIC: 447.46
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_S1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.122
  1. Dissolved oxygen
my_glmS24_O1<- glm(Consumed_individuals~D_oxygen_mgL, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_O1)

AICc(my_glmS24_O1)
## [1] 452.8848
summary(my_glmS24_O1)
## 
## Call:
## glm(formula = Consumed_individuals ~ D_oxygen_mgL, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.39170    0.08033  17.324   <2e-16 ***
## D_oxygen_mgL  0.12467    0.14949   0.834    0.404    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 61.028  on 118  degrees of freedom
## AIC: 452.78
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_O1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.014
  1. Visibility
my_glmS24_V1<- glm(Consumed_individuals~V_Visibility, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_V1)

AICc(my_glmS24_V1)
## [1] 452.1049
summary(my_glmS24_V1)
## 
## Call:
## glm(formula = Consumed_individuals ~ V_Visibility, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.53164    0.08209  18.659   <2e-16 ***
## V_Visibility -0.20381    0.16916  -1.205    0.228    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 60.248  on 118  degrees of freedom
## AIC: 452
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_V1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.030
  1. Community composition
my_glmS24_PC1<- glm(Consumed_individuals~PCo1, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_PC1)

AICc(my_glmS24_PC1)
## [1] 451.8008
summary(my_glmS24_PC1)
## 
## Call:
## glm(formula = Consumed_individuals ~ PCo1, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.50543    0.06171  24.394   <2e-16 ***
## PCo1        -0.16872    0.12776  -1.321    0.187    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.944  on 118  degrees of freedom
## AIC: 451.7
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_PC1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.037
  1. Kelp cover
my_glmS24_K1<- glm(Consumed_individuals~Kelp, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_K1)

AICc(my_glmS24_K1)
## [1] 450.7603
summary(my_glmS24_K1)
## 
## Call:
## glm(formula = Consumed_individuals ~ Kelp, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.37365    0.06285  21.856   <2e-16 ***
## Kelp         0.25829    0.15153   1.705   0.0883 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 58.904  on 118  degrees of freedom
## AIC: 450.66
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_K1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.058
  1. Bare rock cover
my_glmS24_BR1<- glm(Consumed_individuals~Bare_rock, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_BR1)

AICc(my_glmS24_BR1)
## [1] 451.3986
summary(my_glmS24_BR1)
## 
## Call:
## glm(formula = Consumed_individuals ~ Bare_rock, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.51985    0.06556  23.183   <2e-16 ***
## Bare_rock   -0.19684    0.13441  -1.464    0.143    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.542  on 118  degrees of freedom
## AIC: 451.3
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_BR1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.045
  1. Understorey_algae cover
my_glmS24_UA1<- glm(Consumed_individuals~Understorey_algae, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_UA1)

AICc(my_glmS24_UA1)
## [1] 453.5612
summary(my_glmS24_UA1)
## 
## Call:
## glm(formula = Consumed_individuals ~ Understorey_algae, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        1.44308    0.05495  26.261   <2e-16 ***
## Understorey_algae  0.01475    0.12466   0.118    0.906    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 61.705  on 118  degrees of freedom
## AIC: 453.46
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_UA1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.000
  1. Species richness
my_glmS24_SR1<- glm(Consumed_individuals~Mean_SR, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_SR1)

AICc(my_glmS24_SR1)
## [1] 453.2822
summary( my_glmS24_SR1)
## 
## Call:
## glm(formula = Consumed_individuals ~ Mean_SR, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.40627    0.08759  16.054   <2e-16 ***
## Mean_SR      0.08470    0.15635   0.542    0.588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 61.426  on 118  degrees of freedom
## AIC: 453.18
## 
## Number of Fisher Scoring iterations: 4
r2( my_glmS24_SR1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.006
  1. Abundance of invertebrate predators
my_glmS24_I1<- glm(Consumed_individuals~Invertebrate_predators, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_I1)

AICc(my_glmS24_I1)
## [1] 453.4745
summary(my_glmS24_I1)
## 
## Call:
## glm(formula = Consumed_individuals ~ Invertebrate_predators, 
##     family = poisson, data = Squidpops24_stand)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             1.46575    0.07385  19.848   <2e-16 ***
## Invertebrate_predators -0.04313    0.13606  -0.317    0.751    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 61.618  on 118  degrees of freedom
## AIC: 453.37
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_I1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.002
  1. Abundance of predatory fish
my_glmS24_PF1<- glm(Consumed_individuals~Predatory_fish, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_PF1)

AICc(my_glmS24_PF1)
## [1] 451.6715
summary(my_glmS24_PF1)
## 
## Call:
## glm(formula = Consumed_individuals ~ Predatory_fish, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.37803    0.06749  20.418   <2e-16 ***
## Predatory_fish  0.17727    0.12762   1.389    0.165    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.815  on 118  degrees of freedom
## AIC: 451.57
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_PF1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.039
  1. Fish length (cm)
my_glmS24_FL1<- glm(Consumed_individuals~Fish_size, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_FL1)

AICc(my_glmS24_FL1)
## [1] 451.6657
summary(my_glmS24_FL1)
## 
## Call:
## glm(formula = Consumed_individuals ~ Fish_size, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.3041     0.1144  11.400   <2e-16 ***
## Fish_size     0.2283     0.1665   1.371     0.17    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.809  on 118  degrees of freedom
## AIC: 451.56
## 
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_FL1)
## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.039
  • Test interaction between temperature and dissolved oxygen
my_glm<- glm(Consumed_individuals~ Temperature * D_oxygen_mgL, data=Squidpops24_stand, family = poisson)
summary(my_glm)
## 
## Call:
## glm(formula = Consumed_individuals ~ Temperature * D_oxygen_mgL, 
##     family = poisson, data = Squidpops24_stand)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.2170     0.1628   7.474 7.78e-14 ***
## Temperature                0.4416     0.2252   1.961   0.0499 *  
## D_oxygen_mgL              -0.4196     0.4782  -0.877   0.3803    
## Temperature:D_oxygen_mgL   0.3665     0.5390   0.680   0.4965    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 50.734  on 116  degrees of freedom
## AIC: 446.49
## 
## Number of Fisher Scoring iterations: 4
AICc(my_glm)
## [1] 446.8355

Part 6: Herbivory intensity: Effect of site specific variables

First, create datasets to be used in the analysis

Herbivory_analysis <- Herbivory_6 %>%
  dplyr::select(-Rope, -Bait, -Camera, -Hours, -Initial_length_mm, -Final_Length_mm, -Scraped, -Survival_proportion)

1) Standardise

Herb_all_num <- Herbivory_analysis %>%
  dplyr::select(-Site, -Bait2,-Consumed_proportion, -Day_length) 
Herb_all_stand <- decostand(Herb_all_num, method="range")
#Test autocorrelation 
corrplot(cor(Herb_all_stand), method = "number")

Herb_all_stand$Site <- Herbivory_analysis$Site
Herb_all_stand$Consumption<- Herbivory_analysis$Consumption
Herb_all_stand$Consumed_proportion<- Herbivory_analysis$Consumed_proportion

6a) Model selection

  1. Temperature
my_glmL_T1<- glmer(Consumed_proportion~  Temperature + (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_overdispersion(my_glmL_T1)
## # Overdispersion test
## 
##  dispersion ratio = 0.876
##           p-value = 0.584
## No overdispersion detected.
check_model(my_glmL_T1)

AICc(my_glmL_T1)
## [1] 444.0932
summary(my_glmL_T1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ Temperature + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.1    457.2   -219.0    438.1      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.87459 -0.28478 -0.08928  0.12354  3.01044 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.183    2.487   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -0.962      1.740  -0.553     0.58
## Temperature   -2.399      2.515  -0.954     0.34
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temperature -0.883
r2(my_glmL_T1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.671
##      Marginal R2: 0.052
  1. Salinity
my_glmL_S1<- glmer(Consumed_proportion~  Salinity_ppt + (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_S1)

AICc(my_glmL_S1)
## [1] 441.0416
summary(my_glmL_S1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ Salinity_ppt + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    441.0    454.1   -217.5    435.0      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93354 -0.27675 -0.09387  0.13458  2.96049 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 4.546    2.132   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)     1.288      1.846   0.698   0.4853  
## Salinity_ppt   -5.561      2.705  -2.055   0.0398 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Salinty_ppt -0.922
r2(my_glmL_S1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.662
##      Marginal R2: 0.194
  1. Dissolved oxygen
my_glmL_O1<- glmer(Consumed_proportion~  D_oxygen_mgL + (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_O1)

AICc(my_glmL_O1)
## [1] 444.2172
summary(my_glmL_O1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ D_oxygen_mgL + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.2    457.3   -219.1    438.2      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.86805 -0.26809 -0.09095  0.11924  2.97058 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.02     2.454   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    -3.520      1.451  -2.425   0.0153 *
## D_oxygen_mgL    2.400      2.654   0.904   0.3659  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## D_oxygn_mgL -0.831
r2(my_glmL_O1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.664
##      Marginal R2: 0.050
  1. Visibility
my_glmL_V1<- glmer(Consumed_proportion~ V_Visibility + (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_V1)

AICc(my_glmL_V1)
## [1] 444.6448
summary(my_glmL_V1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ V_Visibility + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.6    457.7   -219.3    438.6      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.86293 -0.28102 -0.08826  0.12687  3.02234 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.441    2.538   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    -3.223      1.509  -2.136   0.0327 *
## V_Visibility    1.713      2.851   0.601   0.5480  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## V_Visibilty -0.835
r2(my_glmL_V1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.669
##      Marginal R2: 0.021
  1. Community composition
my_glmL_PC1<- glmer(Consumed_proportion~ PCo1+ (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_PC1)

AICc(my_glmL_PC1)
## [1] 444.2705
summary(my_glmL_PC1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ PCo1 + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.2    457.4   -219.1    438.2      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.88774 -0.28227 -0.08271  0.12286  3.01571 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.24     2.498   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   -3.159      1.175  -2.687  0.00721 **
## PCo1           1.863      2.169   0.859  0.39060   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## PCo1 -0.715
r2(my_glmL_PC1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.670
##      Marginal R2: 0.044
  1. Kelp cover
my_glmL_K1<- glmmTMB(Consumed_proportion~Kelp+ (1|Site), data=Herb_all_stand, family = binomial())              
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
check_model(my_glmL_K1)
## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(my_glmL_K1)
## [1] 451.4504

Linear model failed to converge, quadratic model was chosen for comparisons.

summary(my_glmL_K1)
##  Family: binomial  ( logit )
## Formula:          Consumed_proportion ~ Kelp + (1 | Site)
## Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    451.4    464.5   -222.7    445.4      587 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 5.716    2.391   
## Number of obs: 590, groups:  Site, 12
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -2.1796     1.0412  -2.093   0.0363 *
## Kelp         -0.8131     2.5535  -0.318   0.7502  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(my_glmL_K1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.637
##      Marginal R2: 0.006
  1. Bare rock cover
my_glmL_BR1<- glmer(Consumed_proportion~Bare_rock+ (1|Site), data=Herb_all_stand, family = binomial())         
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_BR1)

AICc(my_glmL_BR1)
## [1] 443.0529
summary(my_glmL_BR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ Bare_rock + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    443.0    456.2   -218.5    437.0      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.88882 -0.26670 -0.09483  0.12657  2.98644 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 5.49     2.343   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   -3.630      1.189  -3.054  0.00226 **
## Bare_rock      3.035      2.136   1.421  0.15545   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Bare_rock -0.758
r2(my_glmL_BR1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.665
##      Marginal R2: 0.107
  1. Understorey algae cover
my_glmL_UA1<- glmer(Consumed_proportion~Understorey_algae+ (1|Site), data=Herb_all_stand, family = binomial())         
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_UA1)

AICc(my_glmL_UA1)
## [1] 437.193
summary(my_glmL_UA1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ Understorey_algae + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    437.2    450.3   -215.6    431.2      587 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8504 -0.2511 -0.0863  0.1267  5.0475 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 3.436    1.854   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        -1.0496     0.7058  -1.487  0.13698   
## Understorey_algae  -6.1317     2.3255  -2.637  0.00837 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Undrstry_lg -0.445
r2(my_glmL_UA1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.714
##      Marginal R2: 0.415
  1. Species richness
my_glmL_SR1<- glmer(Consumed_proportion~Mean_SR+ (1|Site), data=Herb_all_stand, family = binomial())              
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_SR1)

AICc(my_glmL_SR1)
## [1] 444.5391
summary(my_glmL_SR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ Mean_SR + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.5    457.6   -219.2    438.5      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.85322 -0.27247 -0.09101  0.12793  2.97373 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.295    2.509   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -3.364      1.556  -2.162   0.0306 *
## Mean_SR        1.834      2.674   0.686   0.4926  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## Mean_SR -0.849
r2(my_glmL_SR1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.666
##      Marginal R2: 0.027
  1. Abundance of herbivores
my_glmL_H1<- glmmTMB(Consumed_proportion~Herbivores+ (1|Site), data=Herb_all_stand, family = binomial())         
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
check_model(my_glmL_H1)
## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(my_glmL_H1)
## [1] 451.3434
summary(my_glmL_H1)
##  Family: binomial  ( logit )
## Formula:          Consumed_proportion ~ Herbivores + (1 | Site)
## Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    451.3    464.4   -222.7    445.3      587 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 5.732    2.394   
## Number of obs: 590, groups:  Site, 12
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -2.6513     0.9626  -2.754  0.00588 **
## Herbivores    1.1162     2.4550   0.455  0.64934   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(my_glmL_H1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.640
##      Marginal R2: 0.012
  1. Phlorotannins
my_glmL_PL1<- glmer(Consumed_proportion~phloro_mean+ (1|Site), data=Herb_all_stand, family = binomial())              
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_PL1)

AICc(my_glmL_PL1)
## [1] 431.1052
summary(my_glmL_PL1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ phloro_mean + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    431.1    444.2   -212.5    425.1      587 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8432 -0.3082 -0.0483  0.1339 27.0033 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 2.507    1.583   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.08397    0.76929   0.109  0.91309   
## phloro_mean -9.85455    3.39879  -2.899  0.00374 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## phloro_mean -0.689
r2(my_glmL_PL1)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.792
##      Marginal R2: 0.633
  • Test interaction between temperature and dissolved oxygen
my_glm<- glmer(Consumed_individuals~ Temperature*D_oxygen_mgL + (1|Site), data=Crabs2_stand, family = poisson())
summary(my_glm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature * D_oxygen_mgL + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    421.4    435.3   -205.7    411.4      115 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96584 -0.90387  0.06297  0.57760  2.96316 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.2205   0.4696  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)               -0.6119     0.6037  -1.014    0.311
## Temperature                1.0211     0.8251   1.238    0.216
## D_oxygen_mgL              -0.2223     1.7120  -0.130    0.897
## Temperature:D_oxygen_mgL   1.3775     1.9254   0.715    0.474
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr D_xy_L
## Temperature -0.849              
## D_oxygn_mgL -0.735  0.519       
## Tmprtr:D__L  0.730 -0.674 -0.940
AICc(my_glm)
## [1] 421.9314
my_glm<- glmer(Consumed_individuals~ Temperature*Salinity_ppt + (1|Site), data=Crabs24_stand, family = poisson())
## boundary (singular) fit: see help('isSingular')
summary(my_glm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature * Salinity_ppt + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    435.8    449.7   -212.9    425.8      115 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1075 -0.2446  0.1351  0.2772  1.3584 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0        0       
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                0.3196     0.6541   0.489    0.625
## Temperature                0.8991     1.2032   0.747    0.455
## Salinity_ppt               0.8443     1.0544   0.801    0.423
## Temperature:Salinity_ppt  -0.3734     1.7801  -0.210    0.834
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr Slnty_
## Temperature -0.964              
## Salinty_ppt -0.978  0.944       
## Tmprtr:Sln_  0.961 -0.987 -0.972
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
AICc(my_glm)
## [1] 436.2851
my_glm<- glmer(Consumed_individuals~ Temperature*Salinity_ppt + (1|Site), data=Squidpops24_stand, family = poisson())
## boundary (singular) fit: see help('isSingular')
summary(my_glm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature * Salinity_ppt + (1 | Site)
##    Data: Squidpops24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##      448      462     -219      438      115 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8095 -0.1049  0.1330  0.1700  0.9803 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0        0       
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)  
## (Intercept)               0.99054    0.58952   1.680   0.0929 .
## Temperature               0.43957    1.13231   0.388   0.6979  
## Salinity_ppt              0.28427    0.95104   0.299   0.7650  
## Temperature:Salinity_ppt -0.06261    1.66684  -0.038   0.9700  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr Slnty_
## Temperature -0.969              
## Salinty_ppt -0.979  0.948       
## Tmprtr:Sln_  0.966 -0.988 -0.974
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
AICc(my_glm)
## [1] 448.5681

Part 7: Visualisation of best fitted predictors

7a): Crab consumption

Make predictions according to individual models with lowest AICc

  1. Crabs 2 h
p_C2T <- ggpredict(my_glmT1, terms ="Temperature [all]")
p_C2Ox <- ggpredict(my_glmOxC2_1, terms ="D_oxygen_mgL [all]")
  1. Crabs 24 h
p_C24T <- ggpredict(my_glmC24_T1, terms ="Temperature [all]")
p_C24S <- ggpredict(my_glmC24_S1, terms ="Salinity_ppt [all]")
p_C24BR <- ggpredict(my_glmC24_BR1, terms ="Bare_rock [all]")
p_C24K <- ggpredict(my_glmC24_K1, terms ="Kelp [all]")
p_C24FL <- ggpredict(my_glmC24_FL1, terms ="Fish_size [all]")

1) Temperature vs Crabs

ggplot() +
  geom_ribbon(data = p_C2T, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "darkslateblue", alpha = .25) +
  geom_line(data = p_C2T, aes(x = x, y =predicted), color =  "darkslateblue", linewidth = 0.6) +
  geom_jitter(data = Crabs2_stand, aes(x= Temperature, y=Consumed_individuals), color = "darkslateblue", width=0.008, height = 0.1, size = 1, alpha = 0.5)+
  geom_ribbon(data = p_C24T, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_C24T, aes(x = x, y =predicted), color =  "orange1", linewidth = 0.6) +
  geom_jitter(data = Crabs24_stand, aes(x= Temperature, y=Consumed_individuals), color = "orange1", width=0.008, height = 0.1, size = 1, alpha = 0.5)+
  labs(x = "Temperature (°C)", y = "Crabs consumed")  + Theme1 +
  scale_x_continuous(breaks = c(-0.04868, 0.21096, 0.4706, 0.73023, 0.98987), labels = c("10", "12", "14", "16", "18"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) 
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("TempC.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).

2) Salinity vs Crabs

ggplot() +
  geom_ribbon(data = p_C24S, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_C24S, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Crabs24_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
  geom_jitter(data = Crabs2_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
  labs(x = "Salinity (ppt)", y = "Crabs consumed")  + Theme1 +
  scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,8), breaks = seq(0, 6, 1)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("SalC.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).

3) Oxygen vs Crabs

ggplot() +
  geom_ribbon(data = p_C2Ox, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "darkslateblue", alpha = .25) +
  geom_line(data = p_C2Ox, aes(x = x, y = predicted), color =  "darkslateblue", size = 1) +
  geom_jitter(data = Crabs2_stand, aes(x= D_oxygen_mgL, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
  geom_jitter(data = Crabs24_stand, aes(x= D_oxygen_mgL, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
  labs(x = "Dissolved oxygen (mgL)", y = "Crabs consumed")  + Theme1 +
  scale_x_continuous(breaks = c(0.1082296, 0.384854, 0.66147994, 0.935517247), labels = c("4", "6", "8", "10"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,8), breaks = seq(0, 6, 1)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("OxC.jpeg", units="mm", width=90, height=60, dpi=300)
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

4) Bare rock Crabs 24 h

ggplot() +
  geom_ribbon(data = p_C24BR, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_C24BR, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Crabs24_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  geom_jitter(data = Crabs2_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Bare rock cover (%)", y = "Consumed crabs")  + Theme1 +
  scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("BareC.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).

4) Fish length - Crabs 24 h

ggplot() +
  geom_ribbon(data = p_C24FL, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_C24FL, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Crabs24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  geom_jitter(data = Crabs2_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Fish length (cm)", y = "Consumed crabs")  + Theme1 +
  scale_x_continuous(breaks = c(0.011979959
,0.165667696
, 0.319355434
, 0.473043171
, 0.626730908
, 0.780418645
,0.934106383), labels = c("14","16", "18", "20", "22", "24", "26"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("Fish_length.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

7b): Squidpop consumption

  1. Squidpops 2 h

None

  1. Squidpops 24 h
p_S24T <- ggpredict(my_glmS24_T1, terms ="Temperature [all]")
p_S24S <- ggpredict(my_glmS24_S1, terms ="Salinity_ppt [all]")

1) Temperature vs squidpops

ggplot() +
  geom_ribbon(data = p_S24T , aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_S24T , aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Squidpops24_stand, aes(x= Temperature, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
   geom_jitter(data = Squidpops2_stand, aes(x= Temperature, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Temperature (°C)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(-0.04868, 0.21096, 0.4706, 0.73023, 0.98987), labels = c("10", "12", "14", "16", "18"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("TemperatureS.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

2) Salinity vs squidpops

ggplot() +
  geom_ribbon(data = p_S24T, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_S24T, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Squidpops24_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
   geom_jitter(data = Squidpops2_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Salinity (ppt)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("SalS.jpeg", units="mm", width=110, height=80, dpi=300)
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 30 rows containing missing values or values outside the scale range
## (`geom_point()`).
  1. Bare rock vs squidpops
ggplot() +
  geom_jitter(data = Squidpops2_stand, aes(x= Fish_size, y=Consumed_individuals), color = "darkslateblue", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
    geom_jitter(data = Squidpops24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Bare rock cover (%)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("BRS.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
  1. Fish length vs squidpops
ggplot() +
  geom_jitter(data = Squidpops2_stand, aes(x= Fish_size, y=Consumed_individuals), color = "darkslateblue", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
    geom_jitter(data = Squidpops24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Fish Length (cm)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0.011979959
,0.165667696
, 0.319355434
, 0.473043171
, 0.626730908
, 0.780418645
,0.934106383), labels = c("14","16", "18", "20", "22", "24", "26"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("FLS2.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

7c): Seaweed consumption

p_Lphloro <- ggpredict(my_glmL_PL1, terms ="phloro_mean [all]")
p_L24S <- ggpredict(my_glmL_S1, terms ="Salinity_ppt [all]")
p_L24UA <- ggpredict(my_glmL_UA1, terms ="Understorey_algae [all]")
  1. Phlorotannins vs herbivory
ggplot() +
  geom_ribbon(data = p_Lphloro, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_Lphloro, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Herb_all_stand, aes(x= phloro_mean, y=Consumed_proportion), color = "orange1", width=0.02, height = 0.01, size = 1.5, alpha = 0.5)+
  labs(x = "Phlorotannin content (mg"~g^-1*" DW)", y = "Blade consumed proportion")  + Theme1 +
  scale_x_continuous(breaks = c(0.04788192, 0.2840443, 0.5202066, 0.756369, 0.9925314), labels = c("8", "16", "24", "32", "40"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 230 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("PhloroL.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 229 rows containing missing values or values outside the scale range
## (`geom_point()`).
  1. Salinity vs herbivory
ggplot() +
  geom_ribbon(data = p_L24S , aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_L24S, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Herb_all_stand, aes(x= Salinity_ppt, y=Consumed_proportion), color = "orange1", width=0.02, height = 0.01, size = 1.5, alpha = 0.5)+
  labs(x = "Salinity (ppt)", y = "Blade consumed proportion")  + Theme1 +
  scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 268 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("SalL.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 240 rows containing missing values or values outside the scale range
## (`geom_point()`).
  1. UA vs herbivory
ggplot() +
  geom_ribbon(data = p_L24UA, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_L24UA, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Herb_all_stand, aes(x= Understorey_algae, y=Consumed_proportion), color = "orange1", width=0.03, height = 0.03, size = 1.5, alpha = 0.5)+
  labs(x = "UA cover (%)", y = "Predicted consumed proportion ")  + Theme1 +
  scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) + 
  theme(legend.position = c(0.1, 0.85))
## Warning: Removed 269 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("UAL.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 275 rows containing missing values or values outside the scale range
## (`geom_point()`).

Part 8: Validate sea water parameters

Comparison between data obtained in-situ and long-term data. Sea surface temperature (SST) and dissolved oxygen daily average were extracted from E.U. Copernicus Marine Service Information (doi =10.48670/moi-00169, doi =10.48670/moi-00015, respectively). Salinity data were extracted using the NASA Earth Data Search (doi = SMP20-4U7CS). SST and salinity data were collected between 1 December 2018 and 31 March 2019, which matches the period of the study. Dissolved oxygen data were available for the time period 1 November 2021 and 28 February 2025.

8a) Temperature

Temp <- ggplot() +
  geom_jitter(data = temp_data, aes(x= Latitude*-1, y=analysed_sst), width = 0.08, height = 0.03, size = 0.8, alpha = 0.3) +
  geom_boxplot(data = temp_data, aes(x= Latitude*-1, y=analysed_sst, group =  cut_width(Latitude,1)), alpha = 0.2) +
  geom_point(data =Environmental_sum, aes(x= Latitude*-1,y=Environmental_sum$Temperature), size=2, color = "red", pch =15)+
  Theme1 +
labs( y="SST (°C)", x="Latitude (°S)") +
  theme(axis.text.y = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.x =     element_text(size = 11))

Temp

ggsave("TVar.jpeg", units="mm", width=130, height=100, dpi=300)

8b) Salinity

Sal <- ggplot() +
  geom_jitter(data = Sal_data, aes(x= lat*-1, y=Salinity), width = 0.08, height = 0.03, size = 0.8, alpha = 0.6) +
  geom_boxplot(data = Sal_data, aes(x= lat*-1, y=Salinity, group =  cut_width(lat,1)), alpha = 0.2) +
  geom_point(data =Environmental_sum, aes(x= Latitude*-1,y=Environmental_sum$Salinity_ppt), size=2, color = "red", pch =15)+
  Theme1 +
labs( y="Salinity (ppt)", x="Latitude (°S)") +
  theme(axis.text.y = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.x =     element_text(size = 11))

Sal

ggsave("SalVar.jpeg", units="mm", width=130, height=100, dpi=300)

8c) Oxygen

  1. Transform to mg/L
library(respR)
## Warning: package 'respR' was built under R version 4.4.2
## Registered S3 method overwritten by 'respR':
##   method        from
##   print.inspect pryr
## 
## Attaching package: 'respR'
## The following object is masked from 'package:MASS':
## 
##     select
o2_data$conv <- convert_DO(o2_data$o2, # data to convert
                   from = "umol/L",     # oxygen unit to convert from
                   to = "mg/L",    # oxygen unit to convert to
                   t = 15,            # in C
                   S = 34,            # in ppt
                   P = 1.013)         # in bar
  1. Make plot
ox <- ggplot() +
  geom_jitter(data = o2_data, aes(x= Latitude*-1, y=conv), width = 0.08, height = 0.03, size = 0.8, alpha = 0.2) +
  geom_boxplot(data = o2_data, aes(x= Latitude*-1, y=conv, group =  cut_width(Latitude,1)), alpha = 0.2) +
  geom_point(data =Environmental_sum, aes(x= Latitude*-1,y=Environmental_sum$D_oxygen_mgL), size=2, color = "red", pch =15)+
  Theme1 +
labs( y="Dissolved oxygen (mg"~ L^-1*")", x="Latitude (°S)") +
  theme(axis.text.y = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.x =     element_text(size = 11))

ox

ggsave("o2Var.jpeg", units="mm", width=130, height=100, dpi=300)

Part 9: Principal component analysis (PCA) of predictor variables

In this section, a PCA with all the predictor variables was carried out and contracted with latitude.

  1. Get data ready
Variables <- Phlorotannins_sum  %>%
    left_join(Environmental_sum, by = "Site") %>%
  left_join(Group_abundances_mean, by = "Site") %>%
  left_join(Richness_meanpred, by = "Site") %>%
   left_join( pcoa_df1, by = "Site") %>%
   left_join(Size_mean_wide, by = "Site") %>%
   left_join(Benthic_sum_wide_var, by = "Site")

Variables1 <- dplyr::select(Variables, -Site, -Latitude, -Day_length)
  1. Perform PCA
# Perform PCA
pca_result <- prcomp(Variables1, center = TRUE, scale. = TRUE)
# View PCA summary
summary(pca_result)
## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5     PC6     PC7
## Standard deviation     1.9056 1.8232 1.5595 1.2683 1.1188 0.80575 0.67715
## Proportion of Variance 0.2594 0.2374 0.1737 0.1149 0.0894 0.04637 0.03275
## Cumulative Proportion  0.2594 0.4968 0.6705 0.7854 0.8748 0.92121 0.95396
##                            PC8     PC9    PC10    PC11     PC12
## Standard deviation     0.59107 0.42934 0.32003 0.09201 4.97e-15
## Proportion of Variance 0.02495 0.01317 0.00732 0.00060 0.00e+00
## Cumulative Proportion  0.97891 0.99208 0.99940 1.00000 1.00e+00
# Access PCA components
pca_result$rotation  # Loadings (eigenvectors)
##                                PC1         PC2         PC3         PC4
## phloro_mean             0.15312183 -0.38036071  0.16987612 -0.40416794
## D_oxygen_mgL            0.04416129 -0.06947408  0.22900240 -0.29578882
## Salinity_ppt            0.32770211 -0.28725268 -0.30041598  0.01204521
## Temperature             0.33959317 -0.19864160  0.25173653  0.07112323
## V_Visibility           -0.26155336 -0.04655336 -0.46754924 -0.10336321
## Herbivores              0.17945120  0.26495508 -0.32071365 -0.30175664
## Invertebrate_predators  0.23631635  0.35158153  0.08199631 -0.22146569
## Predatory_fish          0.33348022  0.19557271  0.12844257  0.42742255
## Mean_SR                 0.24588446  0.13891112 -0.46658761 -0.10235712
## PCo1                   -0.45336621 -0.12321574  0.20140955  0.11357894
## Fish_size               0.40286210 -0.15837393  0.19130840 -0.19017485
## Kelp                    0.20344908 -0.25734537 -0.07308674  0.49797138
## Bare_rock              -0.02538404  0.41986093  0.34602767 -0.18266636
## Understorey_algae      -0.11785619 -0.43843959 -0.02128394 -0.26379807
##                                PC5          PC6         PC7         PC8
## phloro_mean            -0.22928655 -0.061934912  0.16030316  0.06977640
## D_oxygen_mgL            0.71917764 -0.141970506  0.24870627  0.01063205
## Salinity_ppt           -0.04621847  0.088594787 -0.25384253  0.37791097
## Temperature             0.34233239  0.366339285  0.01239502  0.29709948
## V_Visibility            0.28277131 -0.265296643 -0.32346405 -0.03883120
## Herbivores             -0.18909046 -0.075864110  0.65389541  0.11582022
## Invertebrate_predators  0.07850895 -0.516299347 -0.21984644  0.32372054
## Predatory_fish         -0.04690808 -0.353793051 -0.09874146  0.21766508
## Mean_SR                 0.31766340  0.122121426 -0.05131298 -0.29528597
## PCo1                    0.18079042 -0.201314998  0.07870739  0.07955889
## Fish_size              -0.07309108 -0.101363639 -0.28571931 -0.59602278
## Kelp                    0.08064145 -0.390533765  0.39301185 -0.31521002
## Bare_rock              -0.07793862 -0.008588723 -0.09287449 -0.20915738
## Understorey_algae      -0.18785568 -0.384416932 -0.06340740  0.07047988
##                                PC9         PC10         PC11         PC12
## phloro_mean            -0.19914246 -0.302226749  0.547906843 -0.325191416
## D_oxygen_mgL            0.28917754 -0.203944032 -0.080842398 -0.005764669
## Salinity_ppt            0.29382455 -0.354717840  0.055501156  0.535005519
## Temperature            -0.10134383  0.454952128 -0.004250124 -0.076982913
## V_Visibility            0.07492272 -0.003118584  0.171242875 -0.276713426
## Herbivores              0.13103355 -0.015917972 -0.278008088  0.010513205
## Invertebrate_predators -0.52190028  0.112816781  0.002910929  0.218806222
## Predatory_fish          0.38672379 -0.153625509 -0.023726204 -0.490900382
## Mean_SR                -0.06100320  0.112910544  0.207166319 -0.035055321
## PCo1                   -0.11580691 -0.293793915 -0.143554874  0.226404458
## Fish_size              -0.04289802 -0.178281745 -0.435374536  0.012026704
## Kelp                   -0.11866981  0.127147103  0.290625846  0.278306208
## Bare_rock               0.44193909  0.139799385  0.473690956  0.325001830
## Understorey_algae       0.32715632  0.575421427 -0.145109455  0.028404422
pca_result$x         # Principal components (scores)
##               PC1          PC2         PC3        PC4         PC5        PC6
##  [1,]  1.01934232 -2.465942186 -0.05856192 -1.1805419 -0.74403145 -0.6008792
##  [2,]  1.64364575 -0.247652503  0.33059729  2.8574056  0.41942831 -0.9032240
##  [3,]  0.65673090 -1.475702438 -0.83220277  0.1380973 -0.43291297  1.0595816
##  [4,]  0.10429367  3.053015250 -2.55270043 -1.8166224  0.04174398 -0.5924529
##  [5,]  0.99165334  0.069541318  0.37106187 -0.4478629 -1.35364440  0.5912784
##  [6,] -4.14190953 -0.901118939 -2.06185691  1.0229257 -0.76843221 -0.5500265
##  [7,]  1.57351825  1.793784930  1.24878295 -0.1266223 -0.94011881  0.5183558
##  [8,]  0.52079420  2.322458024  1.45717978  0.2800948  0.28004323 -1.0837141
##  [9,] -0.02668262 -0.009488045 -0.64428859 -0.2340490  2.98842103  0.7665162
## [10,] -0.39575077 -2.747177050  1.78174925 -1.6003581  0.71592698 -0.7627492
## [11,] -3.50506614  1.243557346  2.39160564  0.2027377 -0.10405266  0.8614225
## [12,]  1.55943063 -0.635275706 -1.43136615  0.9047954 -0.10237103  0.6958914
##              PC7         PC8         PC9        PC10        PC11          PC12
##  [1,]  0.2498895 -0.72786447  0.49737574  0.47748733  0.08635639 -4.642120e-15
##  [2,]  0.5267490 -0.03407682 -0.56423696  0.21510138  0.01157954 -5.436623e-15
##  [3,] -0.8074396  0.10588014 -0.19258732  0.36593723 -0.18791464 -5.256212e-15
##  [4,]  0.6350624 -0.05283392 -0.36977732  0.09519209 -0.04917014 -4.992534e-15
##  [5,] -0.6701076 -0.56860844 -0.67721358 -0.40815778  0.11667248 -4.359360e-15
##  [6,] -0.4861734  0.48875365  0.14721247 -0.08360601  0.05770054 -5.280498e-15
##  [7,]  0.1000942  1.32580148  0.29863352  0.19812510  0.07925996 -4.798245e-15
##  [8,] -1.1046837 -0.41762546  0.51534582 -0.17420288 -0.07669713 -4.576201e-15
##  [9,] -0.4462995 -0.02123266  0.04596192  0.08061882  0.10601955 -4.909267e-15
## [10,]  0.3077744  0.66291487 -0.28732577 -0.34383257 -0.06506740 -4.635181e-15
## [11,]  0.8137796 -0.57304152 -0.01998176  0.14372574 -0.03156541 -4.538037e-15
## [12,]  0.8813548 -0.18806684  0.60659323 -0.56638845 -0.04717373 -4.943962e-15
# Scree plot
fviz_eig(pca_result, addlabels = TRUE)

  1. Obtain variables for ggplot
pca_scores <- as.data.frame(pca_result$x)
pca_scores$Site <- Variables$Site
pca_scores$Latitude <- Variables$Latitude
pca_loadings <- as.data.frame(pca_result$rotation)
pca_loadings$Variable <- rownames(pca_loadings)



# Rename variable labels in the pca_loadings data frame
pca_loadings$Variable <- dplyr::recode(pca_loadings$Variable,
  V_Visibility = "Visibility",
  Understorey_algae = "UA cover",
  Bare_rock = "Bare rock cover",
  D_oxygen_mgL = "Dissolved oxygen",
  phloro_mean = "Phlorotannins",
  Kelp = "Kelp cover",
  PCo1 = "PCo1",
  Temperature = "Temperature",
  Fish_size = "Fish length",
  Salinity_ppt = "Salinity",
  Predatory_fish = "Predatory fish",
  Invertebrate_predators = "Invertebrate predators",
  Mean_SR = "Species richness",
  Herbivores = "Herbivores"
)
  1. Visualise
ggplot() +
  # Add PCA scores (samples) with latitude gradient
  geom_point(data = pca_scores, aes(x = PC1, y = PC2, color = Latitude*-1), size = 3) +
  scale_color_viridis_c(option = "plasma", trans = "reverse") +
  # Add PCA loadings (variables as arrows)
  geom_segment(data = pca_loadings,
               aes(x = 0, y = 0, xend = PC1 * 5, yend = PC2 * 5), 
               arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  # Add variable labels
  geom_text(data = pca_loadings, 
            aes(x = PC1 * 5, y = PC2 * 5, label = Variable), 
            color = "black", vjust = -0.5) +
  # Customize plot appearance
  labs(      x = "Principal Component 1",
       y = "Principal Component 2",
       color = "Latitude (°S)") +
  Theme1

ggsave("PCA.jpeg", units="mm", width=130, height=120, dpi=300)
  1. Evaluate statistical effects of latitude
# Load vegan package
library(vegan)

# Create a distance matrix from PCA scores
pca_dist <- dist(pca_scores[, c("PC1", "PC2")])


# Perform PERMANOVA
permanova_result <- adonis2(pca_dist ~ Latitude, data = pca_scores, permutations = 999)
print(permanova_result)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = pca_dist ~ Latitude, data = pca_scores, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)    
## Latitude  1   33.467 0.43742 7.7754  0.001 ***
## Residual 10   43.042 0.56258                  
## Total    11   76.510 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Part 10: Size estimations summary table

A summary table of the size estimations by species and site was created.

Size_summary <- Size_metadata %>%
  group_by(Species) %>%
  summarise (N_individuals = n(), mean_Size_sp = mean(Size_cm), sd = sd(Size_cm), Max_size = max(Size_cm), Min_size = min(Size_cm))
kable(Size_summary)
Species N_individuals mean_Size_sp sd Max_size Min_size
Acanthistius_pictus 3 16.000000 1.0000000 17 15
Anisotremus_scapularis 4 21.250000 2.9860788 25 18
Aplodactylus_punctatus 54 27.944444 6.5196848 40 14
Auchenionchus_critinus 1 12.000000 NA 12 12
Auchenionchus_sp. 14 19.071429 6.0949448 32 11
Auchenionchus_variolosus 3 28.666667 5.6862407 35 24
Bovichtus_chilensis 3 17.666667 9.0184995 27 9
Calliclinus_geniguttatus 19 10.157895 2.3632172 14 4
Cheilodactylus_variegatus 125 22.560000 7.0715240 38 6
Cheilotrema_fasciatum 1 28.000000 NA 28 28
Chromis_crusma 85 13.611765 4.2457096 26 5
Congiopodus_peruvianus 2 13.500000 0.7071068 14 13
Girella_laevifrons 11 26.272727 7.3087743 38 15
Graus_nigra 15 25.066667 4.9923751 37 19
Hemilutjanus_macrophthalmos 2 20.000000 2.8284271 22 18
Hepatus_chilensis 1 6.000000 NA 6 6
Homalaspis_plana 20 9.800000 1.9084301 12 4
Hypsoblennius_sordidus 5 7.800000 1.0954451 9 6
Isacia_conceptionis 7 25.571429 2.5071327 30 22
Labrisomus_philippii 18 28.055556 6.9745101 38 15
Nexilosus_latifrons 4 23.000000 5.7154761 30 17
Paralabrax_humeralis 14 26.571429 8.4826908 44 15
Paralichthys_microps 4 32.500000 23.1012265 65 16
Paraxanthus_barbiger 12 6.250000 1.6583124 9 4
Patagonotothen_brevicauda 10 8.100000 1.3703203 10 6
Pinguipes_chilensis 81 28.370370 8.5197483 51 13
Prolatilus_jugularis 20 24.000000 6.9357958 39 10
Romaleon_setosum 17 8.058824 2.2768012 12 5
Scartichthys_spp. 176 13.670454 5.0927884 34 4
Schroederichthys_chilensis 3 51.666667 2.8867513 55 50
Sebastes_oculatus 12 12.333333 7.1137937 24 4
Semicossyphus_darwini 6 26.500000 6.4420494 36 20
Taliepus_sp. 27 9.444444 1.8045526 12 5
write.csv(Size_summary, "Size_summary2.csv")

Part 11: R Session Information

In this part, information about the operating system and R packages attached during the production of this document can be found

sessionInfo() %>% pander

R version 4.4.1 (2024-06-14 ucrt)

Platform: x86_64-w64-mingw32/x64

locale: LC_COLLATE=English_Australia.utf8, LC_CTYPE=English_Australia.utf8, LC_MONETARY=English_Australia.utf8, LC_NUMERIC=C and LC_TIME=English_Australia.utf8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: respR(v.2.3.3), fitdistrplus(v.1.2-2), survival(v.3.6-4), factoextra(v.1.0.7), lme4(v.1.1-35.4), Matrix(v.1.7-0), glmmTMB(v.1.1.10), mgcv(v.1.9-1), nlme(v.3.1-167), geosphere(v.1.5-18), MuMIn(v.1.48.4), MASS(v.7.3-60.2), corrplot(v.0.95), car(v.3.1-2), carData(v.3.0-5), knitr(v.1.49), performance(v.0.13.0), ggeffects(v.1.7.0), pander(v.0.6.5), patchwork(v.1.3.0), vegan(v.2.6-6.1), lattice(v.0.22-6), permute(v.0.9-7), summarytools(v.1.0.1), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.2), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): rstudioapi(v.0.16.0), DHARMa(v.0.4.6), jsonlite(v.1.9.0), shape(v.1.4.6.1), datawizard(v.1.0.0), magrittr(v.2.0.3), TH.data(v.1.1-2), estimability(v.1.5.1), magick(v.2.8.3), farver(v.2.1.2), nloptr(v.2.1.0), rmarkdown(v.2.29), ragg(v.1.3.2), vctrs(v.0.6.5), minqa(v.1.2.7), base64enc(v.0.1-3), rstatix(v.0.7.2), htmltools(v.0.5.8.1), haven(v.2.5.4), broom(v.1.0.6), marelac(v.2.1.11), sass(v.0.4.9), bslib(v.0.9.0), plyr(v.1.8.9), sandwich(v.3.1-1), emmeans(v.1.10.2), zoo(v.1.8-12), cachem(v.1.1.0), TMB(v.1.9.16), lifecycle(v.1.0.4), pkgconfig(v.2.0.3), sjlabelled(v.1.2.0), R6(v.2.6.1), fastmap(v.1.2.0), rbibutils(v.2.2.16), digest(v.0.6.35), numDeriv(v.2016.8-1.1), colorspace(v.2.1-0), textshaping(v.0.4.0), ggpubr(v.0.6.0), labeling(v.0.4.3), timechange(v.0.3.0), abind(v.1.4-5), compiler(v.4.4.1), withr(v.3.0.2), backports(v.1.5.0), ggsignif(v.0.6.4), gsw(v.1.2-0), roll(v.1.1.7), tools(v.4.4.1), glue(v.1.7.0), grid(v.4.4.1), checkmate(v.2.3.1), cluster(v.2.1.6), reshape2(v.1.4.4), see(v.0.10.0), generics(v.0.1.3), gtable(v.0.3.6), tzdb(v.0.4.0), data.table(v.1.15.4), hms(v.1.1.3), sp(v.2.1-4), xml2(v.1.3.6), ggrepel(v.0.9.6), pillar(v.1.10.1), splines(v.4.4.1), pryr(v.0.1.6), oce(v.1.8-3), tidyselect(v.1.2.1), reformulas(v.0.4.0), stats4(v.4.4.1), xfun(v.0.51), rapportools(v.1.1), matrixStats(v.1.5.0), stringi(v.1.8.4), yaml(v.2.3.8), boot(v.1.3-30), evaluate(v.1.0.3), codetools(v.0.2-20), tcltk(v.4.4.1), seacarb(v.3.3.3), SolveSAPHE(v.2.1.0), cli(v.3.6.2), RcppParallel(v.5.1.10), xtable(v.1.8-4), systemfonts(v.1.1.0), Rdpack(v.2.6), munsell(v.0.5.1), jquerylib(v.0.1.4), Rcpp(v.1.0.12), coda(v.0.19-4.1), parallel(v.4.4.1), bayestestR(v.0.15.1), viridisLite(v.0.4.2), mvtnorm(v.1.2-5), scales(v.1.3.0), insight(v.1.0.2), crayon(v.1.5.3), rlang(v.1.1.4) and multcomp(v.1.4-26)